吳喜之R語言講義(包括各種回歸)

1、吳喜之 免費(fèi)(沒有權(quán)力和銅臭) 資源公開, 可改變代碼(不是黑盒子,也不是吝嗇鬼, 透明是防止“腐敗”的最好方式) 容易學(xué)習(xí)?？删幊桃詫?shí)行復(fù)雜的課題 可擴(kuò)展: 通過數(shù)千個(gè)網(wǎng)上提供的適用于不同領(lǐng)域、不同目的、不同方法的軟件包來實(shí)現(xiàn)你的目標(biāo)。也可以把你的方法貢獻(xiàn)出來 功能強(qiáng)大(繪圖功能, 優(yōu)秀的內(nèi)在幫助系統(tǒng), R社區(qū)的支持，不斷更新，不斷修正) 沒有任何一個(gè)商業(yè)軟件有如此多和如此新的算法 世界應(yīng)用統(tǒng)計(jì)學(xué)家大都把自己的方法首先以R來實(shí)現(xiàn)，并盡量放到R 網(wǎng)站上 一年多，R網(wǎng)站的軟件包數(shù)量增加了兩倍，從近1000個(gè)到近3000多個(gè)。大都都有關(guān)于計(jì)算、演示和輸入輸出方法的函數(shù)和例子數(shù)據(jù) 除非得到巨額資助(

2、或者永遠(yuǎn)使用盜版軟件), 沒有理由在公立學(xué)校教授商業(yè)軟件 絕大多數(shù)美國統(tǒng)計(jì)研究生都會(huì)的語言(Berkeley統(tǒng)計(jì)和應(yīng)用數(shù)學(xué)本科都開設(shè)R語言課) 我的很大一部分?jǐn)?shù)據(jù)分析知識(shí)的來源就是R. 我都能學(xué)會(huì), 并且到處宣傳和普及, 相信你們會(huì)做得更好!點(diǎn)擊點(diǎn)擊CRAN得到一批鏡像網(wǎng)站得到一批鏡像網(wǎng)站下載下載R(/)點(diǎn)擊鏡像網(wǎng)站比如點(diǎn)擊鏡像網(wǎng)站比如Berkeley選擇這個(gè)選擇這個(gè),下載安裝文件下載安裝文件選擇這個(gè)選擇這個(gè),下載軟件包下載軟件包選擇選擇basePackages Packages （每個(gè)都有大量數(shù)據(jù)和可以讀寫修改的（每個(gè)都有大量數(shù)據(jù)和可以讀寫修改的

3、函數(shù)函數(shù)/ /程序）程序）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 （繼續(xù)繼續(xù)）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 （網(wǎng)上）（網(wǎng)上） 網(wǎng)上還有許多所有這些所有這些Packages可以自由下載可以自由下載 Base中的中的package包含常用的函包含常用的函數(shù)和數(shù)據(jù)數(shù)和數(shù)據(jù) 而其他的而其他的packages包含各個(gè)方向包含各個(gè)方向統(tǒng)計(jì)學(xué)家自己發(fā)展的方法和數(shù)據(jù)。統(tǒng)計(jì)學(xué)家自己發(fā)展的方法和數(shù)據(jù)。 希望你是下一個(gè)加盟這些希望你是下一個(gè)加盟這些packages的作者之一。的作者之一。安裝安裝Packages關(guān)機(jī)時(shí)是否保存關(guān)機(jī)時(shí)是否保存? 如果是, 你的運(yùn)算結(jié)果(賦值的變量及函數(shù)等)保存在一個(gè)文件(名字為.RDa

9、ta)中, 下次開機(jī)時(shí)還會(huì)重新載入. 如果你不要?jiǎng)t刪去該文件即可. 其實(shí), 除非是做一個(gè)需要多次才完成的大課題, 一般你都不想保存. 你所用的代碼可以以程序腳本形式(*.R, 注意注意: 一定一定要自己敲入要自己敲入”.R”, 沒有默認(rèn)沒有默認(rèn))保存幾個(gè)有用的函數(shù)函數(shù)：f(x): 名字(變?cè)?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)賦值和運(yùn)算 z = rnorm(1000000,4,0.1) median(z) 賦值賦值: “=”可以用可以用“-”代替代替 xy-w 簡單數(shù)學(xué)運(yùn)算有簡單數(shù)學(xué)運(yùn)算有:+,-,*,/, ,%*%,%(mod) %/%(整數(shù)除法整數(shù)除法)等等等等 常用的數(shù)學(xué)函數(shù)有常用的數(shù)學(xué)函數(shù)有:abs , sign , log , log2, log10 , logb, expm1, log1p(x), sqrt , exp , sin , cos , tan , acos , asin, atan , cosh , sinh， tan

11、h賦值和運(yùn)算 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 , =, , =, &, |, !從高到低的運(yùn)算次序從高到低的運(yùn)算次序一些基本運(yùn)算例子x=1:100(x=1:100)sample(x,20)set.seed(0);sample(1:10

12、,3)#隨機(jī)種子!z=sample(1:200000,10000)z1:10#向量下標(biāo)y=c(1,3,7,3,4,2)zy一些基本運(yùn)算例子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)一些基本運(yùn)算例子pi * 102 #能夠用?”*”來看基本算術(shù)運(yùn)算方法*(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)一些基本運(yùn)算例子class(cars)typeof(cars)names(cars)summary(cars)str(cars)s(cars)class(dist speed)plot(dist speed,cars)一些基本運(yùn)算例子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)一些基本運(yùn)算例子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)一些基本運(yùn)算例子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一些基本運(yùn)算例子#直方圖x - rnorm(200)hist(x, col = light blue)rug(x)#莖葉圖stem(x)#散點(diǎn)圖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)#演示畫圖一些基本運(yùn)算例子#復(fù)數(shù)運(yùn)算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)#行列數(shù)目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的元素的個(gè)數(shù)sum(x,10&x,30|x,1.51,1 #第1中小于.51或者相應(yīng)于第2列中大于0的元素(“或”)x!x,2.51,1#第一列中相應(yīng)于第2列中不小于.51的元素(“非”)apply(x,1,mean);apply(x,2,su

19、m)矩陣/高維數(shù)組#上下三角陣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)得到維數(shù)(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)矩陣/高維數(shù)組/scale#矩陣與向量之間

20、的運(yùn)算 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)#去掉矩陣重復(fù)的行(x=rbind(1:5,runif(5),runif(5),1:5,7:11)x!duplicated(x),unique(x)List#list可以是任何對(duì)象的集合(包括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#上兩行等價(jià)于下面兩行 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)#把a(bǔ)rray變成data.framea=xtabs(FreqSurvived+Sex, w)biplot(corresp(a, nf=2)#應(yīng)用之一 # Start with a data frame.str(mtcar

28、s)x - ftable(mtcarsc(cyl, vs, am, gear)x #為array，其維的次序?yàn)?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)等價(jià)ftable(breakswool+tension,warpbreaks)#as.data.frame(DF - as.data.frame(UCBAdmissions) #等價(jià)于data.frame(UCBAdmissions)xtabs(Freq Admit+ Gender + Dept, DF)#:把方陣變

30、成原來的列聯(lián)表 (a=xtabs( Freq Admit + Gender, data=DF)#如無頻數(shù)(權(quán))，左邊為空寫函數(shù)寫函數(shù)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)#函數(shù)可以不寫return,這時(shí)最后一個(gè)值為return的值.為了輸出多個(gè)值最好使用list#幾個(gè)圖一起: 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)設(shè)置邊緣#缺省值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)關(guān)于畫圖關(guān)于畫圖關(guān)于畫圖關(guān)于畫圖par(mfrow=c(2,2)#準(zhǔn)備畫2x2的4個(gè)圖plot(compression, distance,main= Hookes Law) #只有標(biāo)題plot(co

33、mpression, distance,main= Hookes Law, xlab= x,ylab= y) #標(biāo)題+x,y標(biāo)記identify(compression,distance) #標(biāo)出點(diǎn)號(hào)碼plot(compression, distance,main=Hookes Law) #只有標(biāo)題text(46,120, expression(f=frac(1,2)*k*s)#在指定位寫入文字plot(compression, distance,main=Hookes Law) #只有標(biāo)題的圖text(locator(2), c(I am here!,you are there!) #在點(diǎn)擊

34、的兩個(gè)位置寫入文字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)關(guān)于畫圖關(guān)于畫圖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)關(guān)于畫圖關(guān)于畫圖(符號(hào)顏色大小形狀等)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#調(diào)色板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)關(guān)于畫圖關(guān)于畫圖#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)關(guān)于畫圖關(guān)于畫圖#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關(guān)于畫圖關(guān)于畫圖(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數(shù)據(jù),這里是4個(gè)數(shù)量變量的點(diǎn)圖(最后一個(gè)是分類變量(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數(shù)據(jù),這里只選前7個(gè)數(shù)量變量的點(diǎn)圖#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)關(guān)于畫圖(三維)多窗口操作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()畫圖雜項(xiàng)畫圖雜項(xiàng)#模擬布朗運(yùn)動(dòng)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)#慢動(dòng)作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、散點(diǎn)大小同因變量值成比例x=1:10;y=runif(10)symbols(x,y,circle=y/2,inches=F,bg=x)#數(shù)據(jù)框的每一列都做Q-Q圖table=data.frame(x1=rnorm(100),x2=rnorm(100,1,1)par(ask=TRUE)#waitforchanging等待頁面改變的確認(rèn)results=apply(table,2,qqnorm)par(ask=FALSE)#在一個(gè)圖上添加一個(gè)小圖x=rnorm(100)hist(x)op=par(fig=c(.02,.5,.5,.98),new=TRUE)boxplot(x)#數(shù)學(xué)符號(hào)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統(tǒng)計(jì)模型講義#基礎(chǔ)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()回歸#相關(guān)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#相關(guān)嗎?x=rnorm(3);y=rnorm(3);cor(x,y);cor.test(x,y)$p.valuelibrary(TeachingDemos)put.points.demo()相關(guān)基本原理#基本原理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)#不要截距:等價(jià)于(lm(y-1+x1+x2)summ

54、ary(a);anova(a)names(a)shapiro.test(a$res)qqnorm(a$res);qqline(a$res)#數(shù)學(xué)原理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種處理，用于動(dòng)物實(shí)驗(yàn)，48個(gè)觀測值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