R语言讲义(包括各种回归).ppt

R 语言讲义,吴喜之,免费(没有权力和铜臭)资源公开,可改变代码(不是黑盒子,也不是吝啬鬼,透明是防止“腐败”的最好方式)容易学习。可编程以实行复杂的课题可扩展:通过数千个网上提供的适用于不同领域、不同目的、不同方法的软件包来实现你的目标。也可以把你的方法贡献出来功能强大(绘图功能,优秀的内在帮助系统,R社区的支持，不断更新，不断修正)没有任何一个商业软件有如此多和如此新的算法,世界应用统计学家大都把自己的方法首先以R来实现，并尽量放到R 网站上一年多，R网站的软件包数量增加了两倍，从近1000个到近3000多个。大都都有关于计算、演示和输入输出方法的函数和例子数据除非得到巨额资助(或者永远使用盗版软件),没有理由在公立学校教授商业软件绝大多数美国统计研究生都会的语言(Berkeley统计和应用数学本科都开设R语言课)我的很大一部分数据分析知识的来源就是R.我都能学会,并且到处宣传和普及,相信你们会做得更好!,点击CRAN得到一批镜像网站,下载R(http:/www.r-project.org/),点击镜像网站比如Berkeley,Packages（每个都有大量数据和可以读写修改的函数/程序）,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 Testsforeign 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 Interfacetools Tools for Package Developmentutils The R Utils Package,Packages（继续）,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 effects 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 GeneChipsmaStatistical 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 for Package DevelopmentutilsThe R Utils Package,Packages（网上）,网上还有许多,所有这些Packages可以自由下载,Base中的package包含常用的函数和数据而其他的packages包含各个方向统计学家自己发展的方法和数据。希望你是下一个加盟这些packages的作者之一。,安装Packages,关机时是否保存?,如果是,你的运算结果(赋值的变量及函数等)保存在一个文件(名字为.RData)中,下次开机时还会重新载入.如果你不要则删去该文件即可.其实,除非是做一个需要多次才完成的大课题,一般你都不想保存.你所用的代码可以以程序脚本形式(*.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(1:10)identical(1:10,as.integer(1:10),赋值和运算,z=rnorm(1000000,4,0.1)median(z)赋值:“=”可以用“y-w简单数学运算有:+,-,*,/,%*%,%(mod)%/%(整数除法)等等常用的数学函数有:abs,sign,log,log2,log10,logb,expm1,log1p(x),sqrt,exp,sin,cos,tan,acos,asin,atan,cosh,sinh，tanh,赋值和运算,round,floor,ceilinggamma,lgamma,digamma and trigamma.sum,prod,cumsum,cumprodmax,min,cummax,cummin,pmax,pmin,rangemean,length,var,duplicated,uniqueunion,intersect,setdiff,=,=,&,|,!,从高到低的运算次序,一些基本运算例子,x=1:100(x=1:100)sample(x,20)set.seed(0);sample(1:10,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,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)row.names(cars)class(dist speed)plot(dist speed,cars),一些基本运算例子,head(cars)#cars1:6,tail(cars)ncol(cars);nrow(cars)dim(cars)lm(dist speed,data=cars)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)summaryu(a),一些基本运算例子,x-round(runif(20,0,20),digits=2)summary(x)min(x);max(x)median(x)#medianmean(x)#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 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(x 的最小值=,min(x)#cat(enddocumentn,file=RESULT.tex,append=TRUE)demo(graphics)#演示画图,一些基本运算例子,#复数运算x=2+3i(z-complex(real=rnorm(10),imaginary=rnorm(10)complex(re=rnorm(3),im=rnorm(3)Re(z)Im(z)Mod(z)Arg(z)choose(3,2);factorial(6)#解方程f=function(x)x3-2*x-1uniroot(f,c(0,2)#迭代#如果知道根为极值f=function(x)x2+2*x+1optimize(f,c(-2,2),分布和产生随机数,正态分布:pnorm(1.2,2,1);dnorm(1.2,2,1);qnorm(.7,2,1);rnorm(10,0,1)#rnorm(10)t分布:pt(1.2,1);dt(1.2,2);qt(.7,1);rt(10,1)此外还有指数分布、F分布、“卡方”分布、Beta分布、二项分布、Cauchy分布、Gamma分布、几何分布、超几何分布、对数正态分布、Logistic分布、负二项分布、Poisson分布、均匀分布、Weibull分布、Willcoxon分布等变元可以是向量!,a=factor(letters1:10)a3=w#不行a=as.character(a)a3=wa=factor(a)a,可能遇到的问题,输入输出数据,x=scan()1.5 2.6 3.7 2.1 8.9 12-1.2-4#等价于x=c(1.5,2.6,3.7,2.1,8.9,12,-1.2,-4)setwd(“f:/2010stat”)#或setwd(f:2010stat)(x=rnorm(20)write(x,f:/2010stat/test.txt)y=scan(f:/2010stat/test.txt);yy=iris;y1:5,;str(y)write.table(y,f:/2010stat/test.txt,row.names=F)w=read.table(f:/2010stat/test.txt,header=T)str(w)write.csv(y,f:/2010stat/test.csv)v=read.csv(f:/2010stat/test.csv)str(v)data=read.table(clipboard)write.table(clipboard),序列和向量,z=seq(-1,10,length=100)#z=seq(-1,10,len=100)z=seq(10,-1,-1)#z=10:-1x=rep(1:3,3)x=rep(3:5,1:3)x=rep(c(1,10),c(4,5)w=c(1,3,x,z);w3x=rep(0,10);z=1:3;x+zx*z rev(x)z=c(no cat,has,nine,tails)z1=no catz=1:5z7=8;zz=NULLzc(1,3,5)=1:3;zrnorm(10)c(2,5)z-c(1,3)#去掉第1、3元素z=sample(1:100,10);zwhich(z=max(z)#给出下标,向量矩阵,x=sample(1:100,12);xall(x0);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)#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,10apply(x,2,sum),矩阵/高维数组,#上下三角阵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,#矩阵与向量之间的运算 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.frame,x=matrix(1:6,2,3)z=data.frame(x);zz$X2attributes(z)names(z)=c(TOYOTA,GM,HUNDA)row.names(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,vbindx=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(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)Microbrew 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(beer)/length(beer),Table/categorical data,library(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,lapply,smokes=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(smokes,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 1options(digits=7)#restore the number of digits,array/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(mtcars)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)ftable(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)(DF-as.data.frame(UCBAdmissions)#等价于data.frame(UCBAdmissions)xtabs(Freq Admit+Gender+Dept,DF)#:把方阵变成原来的列联表(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)layout(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,132,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(compression,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!)#在点击的两个位置写入文字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(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=row.names(Animals),adj=1.5)#adj=0 implies left adjusted textplot(log(body),log(brain)identify(log(body),log(brain),row.names(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,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)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,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,sin(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=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=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,green3,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+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 lines: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.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 Topographic 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=peru)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 residual 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)#散点大小同因变量值成比例x=1:10;y=runif(10)symbols(x,y,circle=y/2,inches=F,bg=x)#数据框的每一列都做Q-Q图table=data.frame(x1=rn