R編程筆記一、R數據結構在R中，我們一直都在與對象（object）打交道。存儲類型（mode）實數型（real）：整數（integer）、單精度（single）、雙精度（double）虛數型（complex）：如9+11i字符型（character,string）：如"hello"（單雙引號都行）邏輯型（logical）：TRUE,FALSE（簡寫T,F）函數（function）表達式（expression）2.結構化數據向量（vector）：一列數值或字符矩陣（matrix）：m行×n列（各列之間類型都相同）數據框（dataframe）：類似矩陣，但每一列的數據類型可以不同數組（array）：多維度（不是多變量）列表（list）：有諸多成員雜合在一起，這些成員可以是任意類型，甚至是list本身（及其靈活的數據類型）因子（factor）：分類變量時間序列（ts）：時間序列數據使用變量的時候要特別注意，R對大小寫敏感！產生數據3.1簡單的規則序列>1:10#井號是R的注釋符號[1]12345678910>10:1[1]10987654321>seq(1,10,0.5)#等差數列[1]1.01.52.02.53.03.54.04.55.05.5[11]6.06.57.07.58.08.59.09.510.0>rep(2,5)#重復相同的對象[1]22222>rep(1:3,times=3)#觀察與下例的不同[1]123123123>rep(1:3,each=3)[1]111222333>rep(1:3,1:3)[1]1223333.2產生結構化數據向量很容易用函數c()產生：>x=c(9,1,1)>x[1]911>(x=c('Xie','Yi','Hui'))#為什么用括號?？？？？？？[1]"Xie""Yi""Hui"矩陣用matrix()產生：>matrix(1:10,2)#注意：默認按列排列[,1][,2][,3][,4][,5][1,]13579[2,]246810>matrix(1:10,nrow=2,ncol=5,byrow=T)[,1][,2][,3][,4][,5][1,]12345[2,]678910數據框用data.frame()產生>x=data.frame(1:5,4:8)#把若干個向量合成數據框>xX1.5X4.8114225336447558>x=cbind(x,c('A','B','C','D','E'))#綁上一列字符（按列組合，rbind為按行組合）>xX1.5X4.8c("A","B","C","D","E")114A225B336C447D558E>dimnames(x)#看一下x的行列名[[1]][1]"1""2""3""4""5"[[2]] [1]"X1.5"[2]"X4.8"[3]"c(\"A\",\"B\",\"C\",\"D\",\"E\")">colnames(x)#只看列名[1]"X1.5"[2]"X4.8"[3]"c(\"A\",\"B\",\"C\",\"D\",\"E\")">colnames(x)=c('X1','X2','X3')#改列名>x#這樣看起來就舒服多了X1X2X3114A225B336C447D558E因子用factor()產生列表用list()產生時間序列用ts()產生外部讀入數據剪貼板（clipboard），或者SQL、Access數據庫（RODBC包）等例：D:\x.txt文件"V1""V2""V3""V4""V5"11591317226101418>x=read.table('d:\\x.txt',header=T)>xV1V2V3V4V511591317226101418運算4.1算術運算：+,-,*,/,%%（余數）,%/%（整數商）,^（乘方）>5%%2[1]1>5%/%2[1]24.2邏輯運算：&,|,!（且、或、非）；>,<,>=,<=,=="<-"是賦值符號，x<-9與x=9等價！小于負數加上空格x<-9或者x<(-9)>1==T[1]TRUE>2==T[1]FALSE>0==F[1]TRUE>1:10>5[1]FALSEFALSEFALSEFALSEFALSETRUETRUETRUETRUETRUE5.下標的使用向量、因子、時間序列x[i]；矩陣、數據框x[i,j]x[i,]x[,j]；數組就是根據維度多打幾個逗號而已x[i,j,k,…]；列表要用雙重中括號x[[i]]。特殊的使用：在中括號中用邏輯表達式提取元素或給特定元素賦值，例子如下：>x=matrix(1:20,4,5)>x[,1][,2][,3][,4][,5][1,]1591317[2,]26101418[3,]37111519[4,]48121620>x[x>=2&x<16]#提取數值[1]23456789101112131415>x[x>=2&x<16]=NA#為特定元素賦值>x[,1][,2][,3][,4][,5][1,]1NANANA17[2,]NANANANA18[3,]NANANANA19[4,]NANANA1620數學統計函數6.1基礎：最大值max(),最小值min(),均值mean(),標準差sd(),方差var(),相關系數cor(),求和sum(),積prod(),中位數median(),分位數quantile(),對數log(),指數exp(),排列factorial(),組合choose(),四舍五入round(),向下取整floor(),向上取整ceiling(),總結summary(),6.2矩陣累加cumsum(),秩rank(),排序sort(),倒序rev(),矩陣轉置t(),逆矩陣solve(),特征根eigen(),…6.3統計分布的四大金剛：pnorm(),dnorm(),qnorm(),rnorm(),…(p,d,q,r+分布名稱分別構成：分布函數值、密度函數值、分位數、隨機數，如pf()表示F分布函數值，runif()表示產生均勻分布的隨機數)；？？？？？6.4其他（補充）抽樣sample(),線性模型lm(),廣義線性模型glm(),t檢驗t.test()程序控制流（語句）7.1條件/選擇/循環語句選擇控制if(條件){怎樣怎樣}if(條件){怎樣怎樣}else{又怎樣}ifelse()函數：ifelse(條件,滿足則取我,不滿足則取我)switch()函數：多條件選擇（至少我很少用它）循環for(循環變量in一個序列中){怎樣怎樣}while(某條件滿足則){怎樣怎樣}例：>cond=T>if(cond)winDialog('ok','HelloKitty!')>x=0>for(iin10:1){+x=x+i+}>x[1]55>while(x>=1&winDialog("yesno","I'malittlepig...")!='YES'){+x=x-1+}#這個例子比較整人……二、R統計分析1.線性回歸x=runif(100);y=0.2*x+0.1*rnorm(100)fit=lm(y~x)summary(fit)Call:Call:lm(formula=y~x)Residuals:Min1QMedian3QMax-0.230127-0.0678960.0077060.0540870.249110Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)0.0096390.0196380.4910.625x0.1939760.0327115.9304.55e-08***---Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1Residualstandarderror:0.0979on98degreesoffreedomMultipleR-Squared:0.2641,AdjustedR-squared:0.2566F-statistic:35.17on1and98DF,p-value:4.548e-08plot(x,y);abline(fit)fit=lm(y~x+I(x^2))#高次的回歸，注意I()summary(fit)2.方差分析x=gl(3,50);y=rnorm(150)summary(aov(y~x))DfDfSumSqMeanSqFvaluePr(>F)x21.7430.8720.80530.4489Residuals147159.0921.082boxplot(y~x,notch=T)3.位置參數的檢驗？？？？？？、參數方法t檢驗：t.test()非參數方法：wilcox.test(),kruskal.test(),等等條件檢驗：coin包（ConditionalInference）4.廣義線性模型？？？？？？？三個核心概念：因變量：指數分布族自變量的線性組合連接函數（linkfunction）：用來建立因變量的分布參數與之間的連接用得很廣泛的所謂logistic回歸只不過是GLM的特例：分布族為二項分布、連接函數取logit函數。普通的回歸也是GLM的特例：分布族為正態分布、連接函數為f(x)=xR函數：glm(formula,family=gaussian,data,...)例：counts<-c(18,17,15,20,10,20,25,13,12)outcome<-gl(3,1,9)treatment<-gl(3,3)print(d.AD<-data.frame(treatment,outcome,counts))treatmenttreatmentoutcomecounts111182121731315421205221062320731258321393312summary(glm.D93)5.聚類分析K-Means聚類：kmeans()層次聚類：hclust()專做聚類的包：cluster（里面還有很多種聚類方法）6.列聯表的獨立性檢驗一個有用的函數：table()>table(rpois(100,5),rbinom(100,1,.5))Fisher精確檢驗：fisher.test()卡方檢驗：chisq.test()x=matrix(c(12,5,7,7),nc=2)chisq.test(x)Pearson'sChi-squaredtestwithYates'continuitycorrectiondata:xX-squared=0.6411,df=1,p-value=0.42337.混合效應？？？？？？nlme包lme4包8.分位數回歸quantreg包線性分位數回歸函數rq()等等等等……根據自己的研究需要，可以自己寫一點代碼，或者用別人的包。怎么知道該用什么包呢？——讀書！如：AHandbookofStatisticalAnalysesUsingRMixed-EffectsModelsinSandS-Plus三、作圖用兩天把graphics包的所有函數都看一遍，也就到境界了小例子（入門應該還看不懂，就看個意思吧）火海刀山？y=.1+runif(20,.2,1)xx=c(1,1:20,20)yy=c(0,y,0)plot(xx,yy,type='n',xlab='x',ylab='y')for(iin255:0){yy=c(0,y-(1-i/255)*min(y),0)polygon(xx,yy,col=rgb(1,i/255,0),border=NA)Sys.sleep(.05)}時鐘？x=seq(0,2*pi,.01)par(pty='s')plot(sin(x),cos(x),type='l')n=60for(iin1:n){points(0,0,cex=2)arrows(0,0,cos(2*pi/n*(i-1)),sin(2*pi/n*(i-1)),col='white')arrows(0,0,cos(2*pi/n*i),sin(2*pi/n*i))Sys.sleep(1)arrows(0,0,cos(2*pi/n*(i-1)),sin(2*pi/n*(i-1)),col='white')lines(sin(x),cos(x))}winDialog置信區間：f=function(n=1000,alpha=0.95,rn=50){d=replicate(n,rnorm(rn))m=colMeans(d)#s=apply(d,2,sd)z=qnorm(1-(1-alpha)/2)y0=m-z*1/sqrt(50)y1=m+z*1/sqrt(50)plot(1,xlim=c(0.5,n+0.5),ylim=c(min(y0),max(y1)),type="n",xlab="",ylab="")abline(h=0,lty=2)for(iin1:n){arrows(i,y0[i],i,y1[i],length=0.05,angle=90,code=3,col=ifelse(0>y0[i]&0<y1[i],"blue","red"))points(i,m[i])Sys.sleep(.1)}Sys.sleep(3)}replicate(2,f(100,