吴喜之R语言讲义(包括各种回归)

1、吴喜之 免费(没有权力和铜臭) 资源公开, 可改变代码(不是黑盒子,也不是吝啬鬼, 透明是防止“腐败”的最好方式) 容易学习。可编程以实行复杂的课题 可扩展: 通过数千个网上提供的适用于不同领域、不同目的、不同方法的软件包来实现你的目标。也可以把你的方法贡献出来 功能强大(绘图功能, 优秀的内在帮助系统, R社区的支持，不断更新，不断修正) 没有任何一个商业软件有如此多和如此新的算法 世界应用统计学家大都把自己的方法首先以R来实现，并尽量放到R 网站上 一年多，R网站的软件包数量增加了两倍，从近1000个到近3000多个。大都都有关于计算、演示和输入输出方法的函数和例子数据 除非得到巨额资助(

2、或者永远使用盗版软件), 没有理由在公立学校教授商业软件 绝大多数美国统计研究生都会的语言(Berkeley统计和应用数学本科都开设R语言课) 我的很大一部分数据分析知识的来源就是R. 我都能学会, 并且到处宣传和普及, 相信你们会做得更好!点击点击CRAN得到一批镜像网站得到一批镜像网站下载下载R(/)点击镜像网站比如点击镜像网站比如Berkeley选择这个选择这个,下载安装文件下载安装文件选择这个选择这个,下载软件包下载软件包选择选择basePackages Packages （每个都有大量数据和可以读写修改的（每个都有大量数据和可以读写修改的

3、函数函数/ /程序）程序）base The R Base Packageboot Bootstrap R (S-Plus) Functions (Canty)class Functions for Classificationcluster Cluster Analysis Extended Rousseeuw et al.concord Concordance and reliabilitydatasets The R Datasets PackageexactRankTests Exact Distributions for Rank and Permutation Testsforeig

4、n Read Data Stored by Minitab, S, SAS, SPSS, Stata, Systat, dBase, .graphics The R Graphics PackagegrDevices The R Graphics Devices and Support for Colours and Fontsgrid The Grid Graphics PackageKernSmooth Functions for kernel smoothing for Wand & Jones (1995)lattice Lattice Graphics Interfaceto

5、ols Tools for Package Developmentutils The R Utils PackagePackages （继续继续）MASSMain Package of Venables and Ripleys MASSmethodsFormal Methods and ClassesmgcvGAMs with GCV smoothness estimation and GAMMs by REML/PQLmulttestResampling-based multiple hypothesis testingnlmeLinear and nonlinear mixed effec

6、ts modelsnnetFeed-forward Neural Networks and Multinomial Log-Linear ModelsnortestTests for NormalityoutliersTests for outliersplsPartial Least Squares Regression (PLSR) and Principal Component Regression (PCR)pls.pcrPLS and PCR functionsrpartRecursive PartitioningSAGxStatistical Analysis of the Gen

7、eChipsmaStatistical Microarray AnalysisspatialFunctions for Kriging and Point Pattern AnalysissplinesRegression Spline Functions and ClassesstatsThe R Stats Packagestats4Statistical Functions using S4 ClassessurvivalSurvival analysis, including penalised likelihood.tcltkTcl/Tk InterfacetoolsTools fo

8、r Package DevelopmentutilsThe R Utils PackagePackages （网上）（网上） 网上还有许多所有这些所有这些Packages可以自由下载可以自由下载 Base中的中的package包含常用的函包含常用的函数和数据数和数据 而其他的而其他的packages包含各个方向包含各个方向统计学家自己发展的方法和数据。统计学家自己发展的方法和数据。 希望你是下一个加盟这些希望你是下一个加盟这些packages的作者之一。的作者之一。安装安装Packages关机时是否保存关机时是否保存? 如果是, 你的运算结果(赋值的变量及函数等)保存在一个文件(名字为.RDa

9、ta)中, 下次开机时还会重新载入. 如果你不要则删去该文件即可. 其实, 除非是做一个需要多次才完成的大课题, 一般你都不想保存. 你所用的代码可以以程序脚本形式(*.R, 注意注意: 一定一定要自己敲入要自己敲入”.R”, 没有默认没有默认)保存几个有用的函数函数：f(x): 名字(变元)getwd()setwd(dir = f:/2010stat)#或setwd(f:/2010stat)getwd()x=rnorm(100)ls()?rnorm#或help(rnorm)apropos(“norm“)identical(1:10,1:10)identical(1:10,as.numeric

10、(1:10)identical(1:10,eger(1:10)赋值和运算 z = rnorm(1000000,4,0.1) median(z) 赋值赋值: “=”可以用可以用“-”代替代替 xy-w 简单数学运算有简单数学运算有:+,-,*,/, ,%*%,%(mod) %/%(整数除法整数除法)等等等等 常用的数学函数有常用的数学函数有:abs , sign , log , log2, log10 , logb, expm1, log1p(x), sqrt , exp , sin , cos , tan , acos , asin, atan , cosh , sinh， tan

11、h赋值和运算 round, floor, ceiling gamma , lgamma, digamma and trigamma. sum, prod, cumsum, cumprod max, min, cummax, cummin, pmax, pmin, range mean, length, var, duplicated, unique union, intersect, setdiff , =, , =, &, |, !从高到低的运算次序从高到低的运算次序一些基本运算例子x=1:100(x=1:100)sample(x,20)set.seed(0);sample(1:10

12、,3)#随机种子!z=sample(1:200000,10000)z1:10#向量下标y=c(1,3,7,3,4,2)zy一些基本运算例子z=sample(x,20,rep=T)z(z1=unique(z);length(z1)z=sample(x,100,rep=T)xz=setdiff(x,z)sort(union(xz,z)sort(union(xz,z)=xsetequal(union(xz,z),x)intersect(1:10,7:50)sample(1:100,20,prob=1:100)一些基本运算例子pi * 102 #能够用?”*”来看基本算术运算方法*(pi, (10,

13、2)pi * (1:10)2x - pi * 102xprint(x)(x=pi * 102)pi(1:5)print(x, digits = 12)class(x)typeof(x)一些基本运算例子class(cars)typeof(cars)names(cars)summary(cars)str(cars)s(cars)class(dist speed)plot(dist speed,cars)一些基本运算例子head(cars)#cars1:6,tail(cars)ncol(cars);nrow(cars)dim(cars)lm(dist speed, data = ca

14、rs)cars$qspeed =cut(cars$speed, breaks =quantile(cars$speed), include.lowest = TRUE)names(cars)cars3table(cars3)is.factor(cars$qspeed)plot(dist qspeed, data = cars)(a=lm(dist qspeed, data = cars)summary(a)一些基本运算例子x - round(runif(20,0,20), digits=2) summary(x) min(x);max(x)median(x) # medianmean(x) #

15、 meanvar(x) # variancesd(x) # standard deviation sqrt(var(x)rank(x) # rankorder(x)xorder(x)sort(x)sort(x,decreasing=T)#sort(x,dec=T)sum(x);length(x)round(x)一些基本运算例子fivenum(x) # quantilesquantile(x) # quantiles (different convention)有多种定义quantile(x, c(0,.33,.66,1)mad(x) # normalized mean deviation to

16、 the median (“median average distance“) 可用?mad查看cummax(x)cummin(x) cumprod(x)cor(x,sin(x/20) # correlation一些基本运算例子#直方图x - rnorm(200)hist(x, col = light blue)rug(x)#茎叶图stem(x)#散点图N - 500 x - rnorm(N) y - x + rnorm(N) plot(y x)a=lm(yx)abline(a,col=red)#或者abline(lm(yx),col=red)print(Hello World!)paste(

17、x 的最小值= , min(x)#cat(enddocumentn, file=RESULT.tex, append=TRUE)demo(graphics)#演示画图一些基本运算例子#复数运算x=2+3i(z 0);all(x!=0);any(x0);(1:10)x0diff(x)diff(x,lag=2)x=matrix(1:20,4,5);xx=matrix(1:20,4,5,byrow=T);xt(x)x=matrix(sample(1:100,20),4,5)2*xx+5y=matrix(sample(1:100,20),5,4)x+t(y)(z=x%*%y)z1=solve(z) #

18、 solve(a,b)可以解ax=b方程 z1%*%zround(z1%*%z,14)矩阵nrow(x); ncol(x);dim(x)#行列数目x=matrix(rnorm(24),4,6)xc(2,1),#第2和第1行x,c(1,3) #第1和第3列x2,1 #第2,1元素xx,10,1 #第1列大于0的元素sum(x,10) #第1列大于0的元素的个数sum(x,10&x,30|x,1.51,1 #第1中小于.51或者相应于第2列中大于0的元素(“或”)x!x,2.51,1#第一列中相应于第2列中不小于.51的元素(“非”)apply(x,1,mean);apply(x,2,su

19、m)矩阵/高维数组#上下三角阵x=matrix(rnorm(24),4,6)diag(x)diag(1:5)diag(5)xlower.tri(x)=0#xupper.tri(x)=0;diag(x)=0 x=array(runif(24),c(4,3,2);xis.matrix(x) #可由dim(x)得到维数(4,3,2)is.matrix(x1,)x=array(1:24,c(4,3,2)xc(1,3),x=array(1:24,c(4,3,2)apply(x,1,mean)apply(x,1:2,sum)apply(x,c(1,3),prod)矩阵/高维数组/scale#矩阵与向量之间

20、的运算 x=matrix(1:20,5,4)sweep(x,1,1:5,*)x*1:5 sweep(x,2,1:4,+)(x=matrix(sample(1:100,24),6,4);(x1=scale(x)(x2=scale(x,scale=F); (x3=scale(x,center=F)round(apply(x1,2,mean),14)apply(x1,2,sd)round(apply(x2,2,mean),14);apply(x2,2,sd)round(apply(x3,2,mean),14);apply(x3,2,sd)Data.framex=matrix(1:6,2,3)z=da

21、ta.frame(x);zz$X2attributes(z)names(z)=c(TOYOTA,GM,HUNDA)s(z)=c(2001,2002)Zattach(x)GMdetach(x)GMsapply(z,is.numeric)#apply(z,2,is.numeric)缺失值问题等airqualitycomplete.cases(airquality)#哪一行没有缺失值which(complete.cases(airquality)=F)sum(complete.cases(airquality)na.omit(airquality)#append,cbind,vbin

22、dx=1:10;x12=3(x1=append(x,77,after=5)cbind(1:3,4:6);rbind(1:3,4:6)#去掉矩阵重复的行(x=rbind(1:5,runif(5),runif(5),1:5,7:11)x!duplicated(x),unique(x)List#list可以是任何对象的集合(包括lists)z=list(1:3,Tom=c(1:2, a=list(R,letters1:5),w=hi!)z1;z2z$Tz$T$a2z$T3z$T$wattributes(airquality)#属性!airquality$Ozoneattributes(matrix(

23、1:6,2,3)Categorical dataA survey asks people if they smoke or not. The data is Yes, No, No, Yes, Yesx=c(Yes,No,No,Yes,Yes)table(x);xfactor(x)Barplot:Suppose, a group of 25 people are surveyed as to their beer-drinking preference. The categories were (1) Domestic can, (2) Domestic bottle, (3) Microbr

24、ew and (4) import. The raw data is 3 4 1 1 3 4 3 3 1 3 2 1 2 1 2 3 2 3 1 1 1 1 4 3 1beer = scan() 3 4 1 1 3 4 3 3 1 3 2 1 2 1 2 3 2 3 1 1 1 1 4 3 1barplot(beer) # this isnt correctbarplot(table(beer) # Yes, call with summarized databarplot(table(beer)/length(beer) # divide by n for proportiontable(b

25、eer)/length(beer)Table/categorical datalibrary(MASS)quineattach(quine)table(Age)table(Sex, Age); tab=xtabs( Sex + Age, quine); unclass(tab)tapply(Days, Age, mean)tapply(Days, list(Sex, Age), mean)#apply, sapply, tapply, lapplysmokes = c(Y,N,N,Y,N,Y,Y,Y,N,Y)amount = c(1,2,2,3,3,1,2,1,3,2)(tmp=table(s

26、mokes,amount) # store the tableoptions(digits=3) # only print 3 decimal placesprop.table(tmp,1) # the rows sum to 1 nowprop.table(tmp,2) # the columns sum to 1 now#上两行等价于下面两行 sweep(tmp, 1, margin.table(tmp, 1), /)sweep(tmp, 2, margin.table(tmp, 2), /)prop.table(tmp)#amount # all the numbers sum to 1

27、options(digits=7) # restore the number of digitsarray/matrixtabledata.frame# Start with a contingency table.ftable(Titanic, row.vars = 1:3)ftable(Titanic, row.vars = 1:2)data.frame(Titanic)#把array变成data.framea=xtabs(FreqSurvived+Sex, w)biplot(corresp(a, nf=2)#应用之一 # Start with a data frame.str(mtcar

28、s)x - ftable(mtcarsc(cyl, vs, am, gear)x #为array，其维的次序为(cyl, vs, am, gear) ftable(x, row.vars = c(2, 4)#从x（array）确定表的行变量 # Start with expressions, use table()s dnn to change labelsftable(mtcars$cyl, mtcars$vs, mtcars$am, mtcars$gear, row.vars = c(2, 4), dnn = c(Cylinders, V/S, Transmission, Gears)ft

29、able(vscarb,mtcars)#vs是列，carb是行#或ftable(mtcars$vsmtcars$carb)ftable(carbvs,mtcars) #vs是行，carb是列ftable(mtcars,c(8,11)#和上面ftable(carbvs,mtcars)等价ftable(breakswool+tension,warpbreaks)#as.data.frame(DF - as.data.frame(UCBAdmissions) #等价于data.frame(UCBAdmissions)xtabs(Freq Admit+ Gender + Dept, DF)#:把方阵变

30、成原来的列联表 (a=xtabs( Freq Admit + Gender, data=DF)#如无频数(权)，左边为空写函数写函数ss=function(n=100)z=2;for (i in 2:n)if(any(i%2:(i-1)=0)=F)z=c(z,i);return(z) fix(ss)ss()t1=Sys.time()ss(10000)Sys.time()-t1system.time(ss(10000)#函数可以不写return,这时最后一个值为return的值.为了输出多个值最好使用list#几个图一起: par(mfrow=c(2,4)#par(mfcol=c(2,4)lay

31、out(matrix(c(1,1,1,2,3,4,2,3,4),nr=3,byrow=T)hist(rnorm(100),col=Red,10)hist(rnorm(100),col=Blue,8)hist(rnorm(100),col=Green)hist(rnorm(100),col=Brown)#par(mar = c(bottom, left, top, right)设置边缘#缺省值c(5, 4, 4, 2) + 0.1 (英寸)spring= data.frame(compression=c(41,39,43,53,42,48,47,46),distance=c(120,114,13

32、2,157,122,144,137,141)attach(spring)#(Hookes law: f=.5ks)par(mfcol=c(2,2)plot(distance compression)plot(distance compression,type=l)plot(compression, distance,type=o)plot(compression, distance,type=b)关于画图关于画图关于画图关于画图par(mfrow=c(2,2)#准备画2x2的4个图plot(compression, distance,main= Hookes Law) #只有标题plot(co

33、mpression, distance,main= Hookes Law, xlab= x,ylab= y) #标题+x,y标记identify(compression,distance) #标出点号码plot(compression, distance,main=Hookes Law) #只有标题text(46,120, expression(f=frac(1,2)*k*s)#在指定位写入文字plot(compression, distance,main=Hookes Law) #只有标题的图text(locator(2), c(I am here!,you are there!) #在点击

34、的两个位置写入文字par(mfrow=c(1,1)plot(1:10,sin(1:10),type=l,lty=2,col=4,main=paste(strwrap(The title is too long, and I hate to make it shorter,!#$%&*,width=50),collapse=n)legend(1.2,1.0,Just a sine,lty=2,col=4)关于画图关于画图library(MASS);data(Animals);attach(Animals)par(mfrow=c(2,2)plot(body, brain)plot(sqrt

35、(body), sqrt(brain)plot(body)0.1, (brain)0.1)plot(log(body),log(brain) #或者plot(brainbody,log=xy)par(mfrow=c(1,1)par(cex=0.7,mex=0.7) #character (cex) & margin (mex) expansionplot(log(body),log(brain)text(x=log(body), y=log(brain),labels=s(Animals), adj=1.5)# adj=0 implies left adjusted t

36、extplot(log(body),log(brain)identify(log(body),log(brain),s(Animals)关于画图关于画图(符号颜色大小形状等)plot(1,1,xlim=c(1,7.5),ylim=c(0,5),type=n) # Do not plot pointspoints(1:7,rep(4.5,7),cex=seq(1,4,l=7),col=1:7, pch=0:6)text(1:7,rep(3.5,7),labels=paste(0:6,letters1:7),cex=seq(1,4,l=7), col=1:7)points(1:7,

37、rep(2,7), pch=(0:6)+7) # Plot symbols 7 to 13text(1:7)+0.25, rep(2,7), paste(0:6)+7) # Label with symbol numberpoints(1:7,rep(1,7), pch=(0:6)+14) # Plot symbols 14 to 20text(1:7)+0.25, rep(1,7), paste(0:6)+14) # Labels with symbol number#调色板par(mfrow=c(2,4)palette(); barplot(rnorm(15,10,3),col=1:15)

38、palette(rainbow(15);barplot(rnorm(15,10,3),col=1:15)palette(heat.colors(15);barplot(rnorm(15,10,3),col=1:15)palette(terrain.colors(15);barplot(rnorm(15,10,3),col=1:15)palette(topo.colors(15);barplot(rnorm(15,10,3),col=1:15)palette(cm.colors(15);barplot(rnorm(15,10,3),col=1:15)palette(gray(seq(0,0.9,

39、l=15);barplot(rnorm(15,10,3),col=1:15)palette(grey(seq(0,0.5,l=15);barplot(rnorm(15,10,3),col=1:15)palette(default)par(mfrow=c(1,1)关于画图关于画图#matplotsines=outer(1:20,1:4,function(x, y) sin(x/20*pi*y)matplot(sines, pch = 1:4, type = o, col = rainbow(ncol(sines)#legendx - seq(-pi, pi, len = 65)plot(x, s

40、in(x), type = l, ylim = c(-1.2, 1.8), col = 3, lty = 2)points(x, cos(x), pch = 3, col = 4)lines(x, tan(x), type = b, lty = 1, pch = 4, col = 6)title(legend(., lty = c(2, -1, 1), pch = c(-1,3,4), merge = TRUE), cex.main = 1.1)legend(-1, 1.9, c(sin, cos, tan), col = c(3,4,6), lty = c(2, -1, 1), pch =

41、c(-1, 3, 4), merge = TRUE, bg=gray90)关于画图关于画图#barplot and tablepar(mfrow=c(2,2)tN=table(Ni=rpois(100, lambda=5);tNr=barplot(tN, col=gray)lines(r, tN, type=h, col=red, lwd=2) #- type = h plotting *is* barplotbarplot(tN, space = 1.5, axisnames=FALSE, sub = barplot(., space=0, axisnames = FALSE)#如space

42、=1.5则有稀牙缝barplot(tN, space = 0, axisnames=FALSE, sub = barplot(., space=0, axisnames = FALSE)pie(tN)#pie plotpar(mfrow=c(1,1)#加gridplot (1:3)grid(10, 5 , lwd = 2)dev.set;dev.off;dev.list关于画图关于画图(pairs/三维)#pairs#data(iris)pairs(iris1:4, main = Andersons Iris Data - 3 species, pch = 21, bg = c(red, gr

43、een3, blue)unclass(iris$Species)#iris为150 x5数据,这里是4个数量变量的点图(最后一个是分类变量(iris$Species)#stars#data(mtcars)stars(mtcars, 1:7, key.loc = c(14, 1.5), main = Motor Trend Cars : full stars(),flip.labels=FALSE)#mtcars为32x11数据,这里只选前7个数量变量的点图#perspx - seq(-10, 10, length= 30) y - x f - function(x,y) r - sqrt(x2

44、+y2); 10 * sin(r)/r z - outer(x, y, f) zis.na(z) - 1persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = lightblue)data(volcano)par(mfrow=c(2,2)z - 2 * volcano# Exaggerate the reliefx - 10 * (1:nrow(z) #10 meter spacing(S to N)y - 10 * (1:ncol(z) #10 meter spacing(E to W)# Dont draw the grid lin

45、es : border = NA#par(bg = slategray)persp(x, y, z, theta = 135, phi = 30, col = green3, scale = FALSE,ltheta = -120, shade = 0.75, border = NA, box = FALSE)par(bg= white)#contourrx - range(x - 10*1:nrow(volcano)ry - range(y - 10*1:ncol(volcano)ry - ry + c(-1,1) * (diff(rx) - diff(ry)/2tcol - terrain

46、.colors(12)opar - par(pty = s, bg = lightcyan);par(opar)plot(x = 0, y = 0,type = n, xlim = rx, ylim = ry, xlab = , ylab = )u - par(usr)rect(u1, u3, u2, u4, col = tcol8, border = “red”) #rect画矩形contour(x, y, volcano, col = tcol2, lty = solid, add = TRUE, vfont = c(sans serif, plain) title(A Topograph

47、ic Map of Maunga Whau, font = 4)abline(h = 200*0:4, v = 200*0:4, col = lightgray, lty = 2, lwd = 0.1);par(opar)#image x - 10*(1:nrow(volcano) y - 10*(1:ncol(volcano) image(x, y, volcano, col = terrain.colors(100), axes = FALSE) contour(x, y, volcano, levels = seq(90, 200, by=5), add = TRUE, col = pe

48、ru) axis(1, at = seq(100, 800, by = 100) axis(2, at = seq(100, 600, by = 100) box() title(main = Maunga Whau Volcano, font.main = 4)par(mfrow=c(1,1)关于画图(三维)多窗口操作x11()plot(1:10)x11()plot(rnorm(10)dev.set(dev.prev()abline(0,1)# through the 1:10 pointsdev.set(dev.next()abline(h=0, col=gray)# for the re

49、sidual plotdev.set(dev.prev()dev.off(); dev.off()#- close the two X devices#dev.list()画图杂项画图杂项#模拟布朗运动n=100;x=cumsum(rnorm(100);y=cumsum(rnorm(100);plot(x,y,type=l)x=0;y=0;plot(100,ylim=c(-15,15),xlim=c(-15,15)#慢动作for(i in 1:200)x1=x+rnorm(1);y1=y+rnorm(1);segments(x,y,x1,y1);x=x1;y=y1Sys.sleep(.05)#

50、散点大小同因变量值成比例x=1:10;y=runif(10)symbols(x,y,circle=y/2,inches=F,bg=x)#数据框的每一列都做Q-Q图table=data.frame(x1=rnorm(100),x2=rnorm(100,1,1)par(ask=TRUE)#waitforchanging等待页面改变的确认results=apply(table,2,qqnorm)par(ask=FALSE)#在一个图上添加一个小图x=rnorm(100)hist(x)op=par(fig=c(.02,.5,.5,.98),new=TRUE)boxplot(x)#数学符号x=1:10;

51、plot(x,type=n)text(3,2,expression(paste(Temperature(,degree,C) in 2003)text(4,4,expression(bar(x)=sum(frac(xi,n),i=1,n)text(6,6,expression(hat(beta)=(Xt*X).1*Xt*y)text(8,8,expression(zi=sqrt(xi2+yi2)改变大小写字母x=c(I,am,A,BIG, Cat)tolower(x)toupper(x)R统计模型讲义#基础x=rnorm(20,10)t.test(x,m=9,alt=greater)t.tes

52、t(x1:10,m=9,alt=greater)$p.valuet.test(x,con=.90)$confx=rnorm(30,10);y=rnorm(30,10.1)t.test(x,y,alt=less)library(TeachingDemos)ci.examp()run.ci.examp()vis.boxcox()vis.boxcoxu()回归#相关x=rnorm(20);y=rnorm(20);cor(x,y)cor(x,y,method=kendall);cor(x,y,method=spearman)cor.test(x,y);cor.test(x,y,method=kenda

53、ll);cor.test(x,y,method=spearman)cor.test(x,y,method=kendall)$p.value#相关吗?x=rnorm(3);y=rnorm(3);cor(x,y);cor.test(x,y)$p.valuelibrary(TeachingDemos)put.points.demo()相关基本原理#基本原理set.seed(100)x1=rnorm(100);x2=rnorm(100);eps=rnorm(100)y=5+2*x1-3*x2+epsa=lm(yx1+x2)(lm(y0+x1+x2)#不要截距:等价于(lm(y-1+x1+x2)summ

54、ary(a);anova(a)names(a)shapiro.test(a$res)qqnorm(a$res);qqline(a$res)#数学原理x=cbind(1,x1,x2)dim(x)b=solve(t(x)%*%x)%*%t(x)%*%yba$coe3456780510 xy例1：cross.txt63例1: cross.txtw=read.table(cross.txt,header=T)head(w)plot(yx,w);summary(w)a=lm(yx+z,w)summary(a)anova(a)qqnorm(a$res);qqline(a$res)shapiro.test(

55、a$res)a1=lm(yx*z,w)summary(a1);anova(a1)qqnorm(a1$res);qqline(a1$res)shapiro.test(a1$res)anova(a,a1)library(party)#更简单的方法wt=mob(yx|z,data=w)coef(wt);plot(wt)plot(yx,w);abline(coef(wt)1,col=2);abline(coef(wt)2,col=4)回归方程65Poison ExperimentThe data give the survival times (in 10 hour units) in a 3 x 4

56、 factorial experiment, the factors being (a) three poisons and (b) four treatments. Each combination of the two factors is used for four animals, the allocation to animals being completely randomized.Box, G. E. P., and Cox, D. R. (1964). An analysis of transformations (with Discussion). J. R. Statis

57、t. Soc. B, 26, 211-252. /data/general/poison.html 66例2：poison.txt：3种毒药，4种处理，用于动物实验，48个观测值67Poison1.01.52.02.53.03.54.01.01.52.02.53.01.02.03.04.0Treatment1.01.52.02.53.00.81.01.2Time回归setwd(f:/2010stat)w=read.table(poison.txt,head=T)head(w);tail(w)str(

58、w);summary(w)dim(w)w$Poison=factor(w$Poison)w$Treatment=factor(w$Treatment)pairs(w)#直接回归a=lm(TimePoison*Treatment,w)anova(a)a=lm(Time.,w)anova(a)qqnorm(a$res);qqline(a$res)shapiro.test(a$res)#变换a=lm(1/TimePoison+Treatment,w)anova(a)qqnorm(a$res);qqline(a$res)shapiro.test(a$res)summary(a)逐步回归(step或st

59、epAIC MASS)library(MASS)quine.hi - aov(log(Days + 2.5) .4, quine)quine.nxt - update(quine.hi, . . - Eth:Sex:Age:Lrn)quine.stp - stepAIC(quine.nxt, scope = list(upper = Eth*Sex*Age*Lrn, lower = 1), trace = FALSE)quine.stp$anovacpus1 - cpusattach(cpus)for(v in names(cpus)2:7) cpus1v - cut(cpusv, uniqu

60、e(quantile(cpusv), include.lowest = TRUE)detach()cpus0 - cpus1, 2:8 # excludes names, authors predictionscpus.samp - sample(1:209, 100)cpus.lm - lm(log10(perf) ., data = cpus1cpus.samp,2:8)cpus.lm2 - stepAIC(cpus.lm, trace = FALSE)cpus.lm2$anovaexample(birthwt)birthwt.glm - glm(low ., family = binomial, data = bwt)bi