R实战学习笔记

内置数据结构：元组：元组由不同元素构成，每个元素可以存储不同类型的数据，如字符串、数值、元组。元组是写保护的，元组创建后不能再做任何修改。reload()函数将以前导入过的模块再加载一次。重新加载（reload）包括最初导入模块时应用的分析过程和初始化过程。这样就允许在不退出解释器的情况下重新加载已更改的Python模块。Pow(n,2)对某数求平方pow(4.5-4,2)记（Sqrt()求平方根欧几里得距离sqrt(pow(4.5-4)+pow(2-1,2))Deepcopy（）深拷贝，能拷贝对象内部所有数据和引用。若A深拷贝数据B，则A变，B不变。Copy（）浅拷贝，只是复制数据。若A浅拷贝数据B，则A变，B也变。+表示连接%系统R语言入门R区分大小写字母获取帮助help.start()help(foo)或者？fooexample(foo)data()列出加载包中所有可用示例数据集工作空间工作空间是R用来读取文件和保存结果的默认目录。getwd（）setwd(“D:/rworkspace”)D:/rworkspace必须存或者用dir.create()创建新目录load（）从上一次会话结束地方开始输入source（”myscript.R”）将执行包含在文件myscript.R中的R语句集合输出sink(“myoutput”,append=TRUE,split=TRUE)默认情况下，有一样名字文件，会被覆盖原文件。append将文本追加到文件后，而不是覆盖，split将文件发送至屏幕并且输出文件。不加参数的sink（）仅向屏幕输出结果。图形输出Pdf(“mygraphs.pdf”)输出结果作为输入包包是函数、数据、预编译代码以一种定义完善的格式组成的集合。存包的目录称为库函数.libpath()显示库所在位置library()显示库位置以及有哪些包第一次安装一个包用:install.packages(“gclus”)包（引号是必须的）gclus提供了创建增强型散点图的函数。包只需安装一次，但作者可能更新它update.packages(“包名”)查看已经安装的包installed.packages()列出安装包的版本和依赖关系包的载入：包安装时从cran镜像站点下载它加入库的过程需用到时还得library（gclus）包带有示例，help（package=”包名”）可以输出包的简短描述以及包中函数和数据集名称的列表。选中全部代码再按run否则编译器只会运行最后一行代码数据集对象是可以赋值给变量的任何事物。x<-mtcars[order(mpg),]不加,系统分不清向量矩阵数组数据框data.frame()patientid$age用来选取一个给定数据框中的某个特定变量attach(mtcars)//可以将数据框添加到R的搜索路径，plot(disp,mpg)//此时若有其他名称相同对象，比如也叫disp则原始对象取得优先权。若长度不匹配则报错，要注意detach（mtcars）否则得plot(mtcars$disp,mtcars$mpg)with（{只写半边括号，r中命令就能按enter后自动换行命令<<-按上下键可以使用历史命令列表（向量、数组、数据框、列表）因子：因子(factors)提供了一种处理分类数据的更简介的方式。levels=c(“poor”,””,”improved”,”excellent”)可以通过level将类别变量排序列联表table(patientdata$diabete,patientdata$status)str(patientdata)显示数据框patientdata的结构，如哪些列，多少个变量及观测，具体值是。Summary（patient）显示统计概要#注释无多行注释图形初阶如何创建多个图形并随时查看每一个呢方法一dev.new()创建图语句dev.new()创建图语句dev.off()将输出返回到终端修改图形参数>dose<-c(20,30,40,45,60)>drugA<-c(16,20,27,40,60)>drugB<-c(15,18,25,31,40)>plot(dose,drugA,type="b")复制当前图形参数opar<-par(no.readonly=TRUE)par(lty=2,pch=17)#设置图形参数，此时应用于所有之后创建图形plot(dose,drugA,type="b")#查看修改参数后图片效果>par(opar)>plot(dose,drugA,type="b")#图片又变回原来参数了>plot(dose,drugA,type="b",lty=2,pch=17)#仅应用于该图像绘制点的符号R中查看已经安装的包路径和包名library()pchAxis坐标轴font字体类型：1常规2粗体3斜体4粗斜体图形前景色及背景色图形尺寸：代码应该提前，不要等到图已经做出来了再写该代码Par(pin=c(4,3),m,mar=c(1,.5,1,.2))pin长宽Mar指边界大小par(mar=(5,4,4,8)+0.1)下开始，逆时针转注意：某些高级绘图函数已经包含了默认的标题和标签。你可以通过在plot（）语句或者单独的par（）中添加ann=FALSE来移除它们announcement通告坐标轴yaxt="n"xaxt=”n”将分别禁用x轴、y轴(会留下框架，只是禁用了刻度)参考线Abline图例Legend文本标注Textmtexttext(wt,mpg,s(mtcars),cex=.6,pos=4,col="red")直接点的位置，再加上标注向量text(3,3,”exzample”)单独点加标注图形组合Par(mfrow=c(2,2))2行2列按行填充4个图layout(matrix(c(1,1,2,3),2,2,byrow=TRUE),widths=c(3,1),heights=c(1,2))c()为矩阵元素2,2为2行2列layout(matrix(c(1,2,3,0,2,3,0,0,3),nr=3))matrix有9个元素，具有这样的形式：

[,1][,2][,3][1,]

1

0

0[2,]

2

2

0[3,]

3

3

3把这个矩阵传入layout函数，我们就能得到这样的outputdevicepar(fig=c(0,0.8,0,0.8))>plot(wt,mpg,xlab="mile",ylab="weight")par(fig=c(0,0.8,0.55,1),new=TRUE)boxplot(wt,horizontal=TRUE,axes=FALSE)opar<-par(no.readonly=TRUE)par(fig=c(0,0.8,0,0.8))plot(wt,mpg,xlab="mile",ylab="weight")par(fig=c(0,0.8,0.55,1),new=TRUE)boxplot(wt,horizontal=TRUE,axes=FALSE)六、基本图形：条形图barplot()一个数值堆砌条形图分组条形图均值条形图可以使用数据整合函数并将结果传递给barplot()条形图的微调棘状图棘状图对堆砌条形图进行了重缩放，这样每个条形的高度均为1，每一段高度表示比例。library(vcd)attach(Arthritis)counts<-table(Treatment,Improved)spine(counts,main="SpinogramExample")#spine属于vcd包中detach(Arthritis)饼图opar<-par(no.readonly=TRUE)par(mfrow=c(2,2))x<-c(10,12.9,4,16,8)lbls<-c("US","France","Germany","Australia","UK")pie(x,lbls,main="简单饼图")lbls<-paste(lbls,"",round(x/sum(x)*100,2),"%",sep="")#paste()实现字符串连接pie(x,lbls,main="使用百分比标签的饼图",col=rainbow(length(lbls)))install.packages("plotrix")library(plotrix)pie3D(x,labels=lbls,explode=0.1,main="3D饼图")mytable<-table(state.region)lbls<-paste(names(mytable),"\n",mytable,sep="")pie(mytable,labels=lbls,main="列表数据饼图\n(一个简单的)")扇形图比较各个部分大小par(mfrow=c(1,1))

lbls

<-

c("US",

"France",

"Germany",

"Australia",

"UK")#lblslabellist

fan.plot(x,

labels=lbls,

main="扇形图")

直方图histogrambreaks=12,表示分组数hist（x）hist(mtcars$mpg,breaks=12,col="red",xlab="milespergallon",main="coloredhistogramwith12bins")分成12组#根据概率密度画图#增加轴须图hist(mpg,freq=FALSE,#fre=FALSE表示根据概率密度而不是频数绘制图形breaks=12,##

breaks用于控制组的数量col="red",main="直方图、轴须图、概率密度曲线",xlab="MilesPerGallon",ylab="概率密度")rug(jitter(mpg))增加轴须图#rug(jitter(mpg,amount=0.01))如果数据中有许多结（数据中相同的值），可以使用该代码将轴须图的数据打散lines(density(mpg),col="blue",lwd=2)#概率密度曲线核密度曲线轴须图：是实际数据值的一种一维呈现方式x<-mpgh<-hist(x,breaks=12,col="red",main="添加正态密度曲线和外框的直方图",xlab="MilesPerGallon",ylab="频数")xfit<-seq(min(x),max(x),length=40)#生成一个序列，最小最大值为···长度为40个yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))#生成一个yfit<-yfit*diff(h$mids[1:2])*length(x)不理解lines(xfit,yfit,col="blue",lwd=2)box()#外框detach(mtcars)par(opar)核密度图核密度估计是用于估计随机变量概率密度函数的一种非参数方法。Plot(density(x))可比较核密度图箱线图attach(mtcars)boxplot(mpg,main="Boxplot",ylab="milespergallon")并列箱线图跨组比较boxplot(formula,data=data.frame)formula为一个公式，Y~A为类别型变量A每个值并列生成数值型变量y的箱线图Y~A*B为类别型变量A和B所有水平的两两组合生成数值型变量y的箱线图varWidth=TRUE箱线图宽度及样本大小的平方根成正比。boxplot(mpg~cyl,data=mtcars,main="Boxplot",xlab="numberofcylinder",ylab="milespergallon")boxplot(mpg~cyl,data=mtcars,notch=TRUE,varwidth=TRUE,main="Boxplot",xlab="numberofcylinder",ylab="milespergallon")notch=TRUE得到含凹槽的箱线图，若两个箱的凹槽互不重叠则表明他们的中位数显著差异。两个交叉因子的箱线图cyl.f气缸468am.f,自动及标准型号cyl.f<-factor(cyl,levels=c(4,6,8),labels=c("4clinder","6cylinder","8cylinder"))>am.f<-factor(am,levels=c(0,1),labels=c("auto","standard"))>boxplot(mpg~cyl.f*am.f,data=mtcars,varwidth=TRUE,col=c("gold","darkgreen"),main="Boxplot",xlab="numberofcylinder",ylab="milespergallon")am.f*cyl.f这个顺序影响排序小提起图是箱线图和核密度图的合体‘散点图Plot（）点图dotchart()attach(mtcars)x<-mtcars[order(mpg),]x$cyl<-factor(x$cyl)hex$color[x$cyl==4]<-"red"x$color[x$cyl==6]<-"blue"x$color[x$cyl==8]<-"darkgreen"dotchart(x$mpg,labels=s(x),cex=.7,groups=x$cyl,gcolor="black",color=x$color,pch=19,main="gasmileageforcarmodels\ngroupedbycylinder",xlab="milespergallon")groups=x$cyl,分组x$color字符型向量color被添加到数据框中创建分组因子factor(x=character(),levels,labels=levels,exclude=NA,ordered=is.ordered(x),nmax=NA)七、描述性统计分析总体描述library(psych)attach(mtcars)vars<-c("mpg","hp","wt")describe(mtcars[vars])detach(mtcars)结果分组描述性统计量library(reshape)dstats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x)))dfm<-melt(mtcars,measure.vars=c("mpg","hp","wt"),id.vars=c("am","cyl"))cast(dfm,am+cyl+variable~.,dstats)#黄色为类别型分组变量，右边为特征统计量结果列联表考察交叉类型的频数、比例、以及交叉类型之间的相关性一维>mytable<-with(Arthritis,table(Improved))>mytableImprovedNoneSomeMarked421428二维install.packages("gtools")install.packages("gdata")install.packages("gmodels")library(gmodels)with(Arthritis,CrossTable(Treatment,Improved))Help(crosstable)可以计算小数位数、卡方等独立性检验、计算期望皮尔逊残差等。三维其他prop.table()计算比例Margin.table()计算频数独立性、相关性考察数据分析一文件八、回归0LS（ordinaryleastsquare）普通最小二乘法的回归简单线性回归残差标准误可以理解为:模型用身高预测体重的平均误差。R平方：=0.991表明模型可以解释99.1%的方差，它也是实际及预测值之间的相关系数。F统计量：等同于身高回归系数的t检验分析所以可以用一个弯曲曲线来提高预测的精度。多项式回归P<0.001水平下，回归系数均非常显著。模型方差解释率达到99.9%。二次项显著性（t=13.891,p<0.001）表明二次项提高了模型的拟合度。（这样就可以拒绝原假设二次项系数为0的假设）拟合曲线Weight=261.88-7.35heght+0.083height^2多元线性回归预测变量不止一个。例如考察人口、文盲率、收入、霜冻综合对谋杀率的影响。案例：检测变量关系，相关系数矩阵>states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])>cor(states)MurderPopulationIlliteracyIncomeFrostMurder1.00000000.34364280.7029752-0.2300776-0.5388834Population0.34364281.00000000.10762240.2082276-0.3321525Illiteracy0.70297520.10762241.0000000-0.4370752-0.6719470Income-0.23007760.2082276-0.43707521.00000000.2262822Frost-0.5388834-0.3321525-0.67194700.22628221.0000000做散点图矩阵>library(car)scatterplotMatrix(states,spread=FALSE,lty.smooth=2,main="scatterplotmatrix")从图中可以看出谋杀率是双峰曲线，谋杀率随着人口文盲率增加而增加，随着收入水平和结霜天数增加而下降。越冷地方文盲率越低，收入水平越高。多元线性回归参数查看DFdegreeoffreedom自由度有交互项的多元线性回归许多有趣的研究都会涉及交互的预测变量。见书R实战170页用图形展示上图蓝色文字含义回归诊断对模型的信赖程度依赖于它在多大程度上满足OLS模型的统计假设。OLS回归的统计假设正态性:QQ图独立性：因变量是否独立。例如Weight~height。一位女性体重会影响另一位女性体重？如果来自同一个家庭有可能线性同方差性模型在多大程度上满足统计假设的任何信息。标准方法#8.3.1回归诊断标准方法fit<-lm(weight~height,data=women)par(mfrow=c(2,2))plot(fit)上图没有显示因变量是否独立的图改进的方法进行回归诊断Car包提供了大量函数，增强了对拟合回归模型评价的能力。正态性方法一：检验正态性par(mfrow=c(1,1))library(car)states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)qqPlot(fit,labels=s(states),id.method="identify",simulate=TRUE,main="Q-QPlot")#id.method="identify'能够交互式绘图，鼠标单击图形内的点，将会标注函数中label选项的设定值。#simulate=TRUE95%的置信区间将会用参数自助法关注下离群点方法二：残差图residplot<-function(fit,nbreaks=10){z<-rstudent(fit)hist(z,breaks=nbreaks,freq=FALSE,xlab="StudentizedResidual",main="DistributionofErrors")rug(jitter(z),col="brown")curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=2)lines(density(z)$x,density(z)$y,不好理解col="red",lwd=2,lty=2)legend("topright",legend=c("NormalCurve","KernelDensityCurve"),lty=1:2,col=c("blue","red"),cex=.7)residplot(fit)方法三：标准方法中的QQ图，看是否在一条45度直线上。误差的独立性因变量值（或残差）是否相互独立，最好依据收集数据方式的先验知识。比如：时间序列呈现自相关性，相隔时间越近的观测相关性大于相隔越远的观测。Car包中Durbin-Watson函数作Durbin-Watson检验

线性通过成分残差图即偏残差图，可以看看因变量及自变量之间是否线性相关。library(car)crPlots(fit)绿色线为平滑拟合曲线。同方差性Car包中提供了ncvTest()和spreadLevelPlot（）函数判断方差是否恒定。ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性，即误差方差不恒定。spreadLevelPlot（）创建一个添加了最佳拟合曲线的散点图。展示标准化残差绝对值及拟合值的关系。spreadLevelPlot(fit)线性模型假设的综合验证library(gvlma)gvmodel<-gvlma(fit)summary(gvmodel)多重共线性回归系数测量的是当其他预测变量不变时，某个预测变量对响应变量的影响。比如：假定年龄不变，测量握力及年龄的关系？。这种问题称为多重共线性。在回归模型牵涉到多个自变量的时候，自变量之间可能会相互关联，即他们之间存在有多重共线性。它会导致模型参数的置信区间过大，使得单个系数解释起来很困难。方差膨胀因子：VIF（varianceinflationfactor）VIF的平方根表示变量回归参数的置信区间能膨胀为及模型无关的预测变量的程度。异常观测值全面回归分析需要包含对异常值的分析。离群点模型预测效果不佳的观测点，残差比较大的点。判断：1：Q—Q图中落在置信区间以外的点被认为是离群点2：标准化残差值大于2或者小于-2可能是离群点Car包中outlierTest()函数可以求得最大标准化残差绝对值Bonferroni调整后的P值。P=0.048该函数只是根据单个最大残差值得显著性判断是否有离群点。若不显著则没有离群点，若显著则需要删除离群点，再检验是否有其他离群点存在高杠杆值点？及响应变量值没有关系，由许多异常的预测变量值组合起来的。通过帽子统计量判断。hat.plot<-function(fit){p<-length(coefficients(fit))包括截距项若四个预测变量则为5n<-length(fitted(fit))观测数量，比如有50个样本值plot(hatvalues(fit),main="IndexPlotofHatValues")abline(h=c(2,3)*p/n,col="red",lty=2)#水平线为帽子均值2倍和3倍的位置。identify(1:n,hatvalues(fit),names(hatvalues(fit)))hat.plot(fit)p/n即为帽子均值。lm(Murder~Population+Illiteracy+Income+Frost,data=states)谋杀率大小由4个预测变量及一个截距（常数项）分别贡献。强影响点对模型参数估计值影响有些比例失衡的点。比如若移除模型的一个观测点时，模型会发生巨大的变化，则需要检测数据中是否存在强影响点。方法：cook距离（D统计量）、CookD值大于4/(n-k-1)则表明她是强影响点。N为样本量大小，k为预测变量的数目。cutoff

<-

4/(nrow(states)-length(fit$coefficients)-2)

plot(fit,

which=4,

cook.levels=cutoff)

abline(h=cutoff,

lty=2,

col="red")方法2：变量添加图：上图可以鉴别强影响点，但是并不提供这些点如何影响模型的参数的信息。library(car)avPlots(fit,ask=FALSE,onepage=TRUE,id.method="identify")主要看红线2边是否只有一边有点，另外一边没有点。将离群点、杠杆值和强影响点整合在一张图画中：influencePlot(fit,id.method="identify",main="InfluencePlot",sub="CirclesizeisproportionaltoCook'sdistance")id.method="identify"为交互作图，点击感兴趣的点，即显示信息。对离群点的改进措施如何处理违背回归假设的问题。删除观测点如果推断离群点是数据记录错误等，可以选择将最大离群点、强影响点删除。删除后重新拟合，若仍然有离群点及强影响点，则重复上面步骤，直至获得满意结果。若不了解产生离群点原因，则需谨慎对待，不要盲目删除。变量变换当模型不符合正态性、线性、或者同方差性假设时，可以将一个或者多个变量变换。变换多用Yλ代替Yλ常用取值-2

-1

-0.5

0

0.5

1

2若Y是比例数，通常使用logit变换[ln(Y/1-Y)]当违反了正态性假设常常对响应变量Y尝试某种变换。Box-cox正态变换当违反了线性假设时，对预测变量X进行变换常常会比较有用。当违反方差稳定性，可以变换响应变量。具体见上文spreadLevelPlot()谨慎对待变量变换：如果变换了变量，你的解释必须机遇变换后的变量，而不是初始变量。如果变换得没有意义，则应该避免这样做。比如：收入的对数变换好理解。但是自杀意念的频率域抑郁症程度的立方根之间的关系在现实中很难解释出有意义的含义。添加或者删除变量是处理多重共线性的非常重要的方法。仅仅做预测则多重共线性并不构成问题，但是如果要解释每个预测变量，则必须解决该问题。常用方法：删除某个存在多重共线性的变量（vif>2岭回归—多元回归的变体使用其他回归模型如果存在离群点或者强影响点可以使用稳健回归模型替代OLS回归模型。如果违背正态性假设可以用非参数回归模型。如果存在显著非线性，尝试用非线性回归模型。如果违背了误差独立性准则，尝试用专门研究误差结构的模型，如：时间序列模型，或者多层次回归模型，或者广义线性模型。选择“最佳”回归模型预测精度及模型简洁度的调和问题。最佳是一种权衡，没有绝对最佳模型。模型比较比较两个嵌套模型的拟合优度。嵌套式指该模型的一些项，完全包含在另一个模型中。在States多元回归模型中可以发现Income和frost回归系数不显著，检验不含这两个变量的模型预测效果是否更好。方法一：anova()函数方法二：AIC赤池信息准则ANOVA需要嵌套模型，而AIC方法不需要变量选择从大量变量中选择最终的预测变量逐步回归模型会一次添加或者删除一个变量，直到达到某个判停准则为止。向前逐步回归：每次添加一个变量向后逐步回归：每次删除一个变量向前向后逐步回归（简称逐步回归）：变量每次进入一个，对贡献不大变量再删除。逐步回归也许找到个好模型，但是不能保证模型就是最佳模型，因为不是每一个模型都被评价了，故有了全子集回归法。Intercept截距项疑问：AIC怎么计算的？全子集回归所有可能的模型都会被检验。全子集回归可用leaps包中regsubsets（）函数实现，通过R平方，调整R平方或者mallowscp统计量等准则来选择最佳模型。R平方含义是预测变量解释响应变量的程度。一般R平方会随着预测变量的增加而增加，但是当及样本量相比，预测变量数目很大时，容易导致过拟合。R平方很可能会丢失数据的偶然变异信息，而调整的R平方提供了更为真实的R平方估计。install.packages("leaps")

library(leaps)

leaps

<-

regsubsets(Murder

~

Population

+

Illiteracy

+

Income

+

Frost,

data=states,

nbest=4)

plot(leaps,

scale="adjr2")mallowscp统计量也是作为逐步回归的判停准则。好的回归模型它的cp统计量非常接近模型的参数数目（包含截距项）library(car)

subsets(leaps,

statistic="cp",

main="Cp

Plot

for

All

Subsets

Regression")

abline(1,

1,

lty=2,

col="red")注意：该图也是交互作图，有时运行没有停止运行，点击下图看是否为大部分情况全子集回归优于逐步回归，但是数据量大时全子集回归会很慢。拟合效果佳但是没有意义的模型没有意义。主题背景知识的理解才能指导获得理想的模型。深层分析模型泛化能力和相对重要性交叉验证一定比例数据作训练集，一定比例数据作检验模型样本。K重交叉验证：样本被分为K个子样本，轮流将K-1个子样本组合为训练集，另外一个子样本作为保留集合。得到K个预测方程，记录K个保留样本的预测表现结果，然后求其平均值。#

代码清单8-15

R平方的k重交叉验证

shrinkage

<-

function(fit,

k=10){

require(bootstrap)

theta.fit

<-

function(x,y){lsfit(x,y)}

theta.predict

<-

function(fit,x){cbind(1,x)%*%fit$coef}

x

<-

fit$model[,2:ncol(fit$model)]

y

<-

fit$model[,1]

results

<-

crossval(x,

y,

theta.fit,

theta.predict,

ngroup=k)

r2

<-

cor(y,

fit$fitted.values)^2

r2cv

<-

cor(y,

results$cv.fit)^2

cat("Original

R-square

=",

r2,

"\n")

cat(k,

"Fold

Cross-Validated

R-square

=",

r2cv,

"\n")

cat("Change

=",

r2-r2cv,

"\n")

}#K重交叉验证函数

fit1

<-

lm(Murder

~

Population

+

Income

+

Illiteracy

+

Frost,

data=states)

shrinkage(fit1)

fit2

<-

lm(Murder

~

Population

+

Illiteracy

,

data=states)

shrinkage(fit2)相对重要性哪些变量对预测最为重要。根据相对重要性对预测变量排序。方法一：预测变量若不相关，则根据预测变量及响应变量的相关系数即可确认重要性排序。但大部分情况下预测变量之间是有一定相关性。方法二标准化的回归系数，表示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化。方法三：相对权重，测量对R平方的贡献附录：#

代码清单8-16

relweights()函数，计算预测变量的相对权重

relweights

<-

function(fit,

...){

R

<-

cor(fit$model)

nvar

<-

ncol(R)

rxx

<-

R[2:nvar,

2:nvar]

rxy

<-

R[2:nvar,

1]

svd

<-

eigen(rxx)

evec

<-

svd$vectors

ev

<-

svd$values

delta

<-

diag(sqrt(ev))

lambda

<-

evec

%*%

delta

%*%

t(evec)

lambdasq

<-

lambda^2

beta

<-

solve(lambda)

%*%

rxy

rsquare

<-

colSums(beta^2)

rawwgt

<-

lambdasq

%*%

beta^2

import

<-

(rawwgt/rsquare)

*

100

lbls

<-

names(fit$model[2:nvar])

rownames(import)

<-

lbls

colnames(import)

<-

"Weights"

barplot(t(import),

names.arg=lbls,

ylab="%

of

R-Square",

xlab="Predictor

Variables",

main="Relative

Importance

of

Predictor

Variables",

sub=paste("R-Square=",

round(rsquare,

digits=3)),

...)

return(import)

}九、方差分析基本术语ANOVA：(AnalysisofVariance方差分析:关注组别差异的分析。组间因子：每个患者仅被分配到一个组别中。组内因子：如植物类型和二氧化碳浓度。每个植物只能划分到一个植物类型A或者B中。但是所A和B中所有植物均可以在相同浓度下测试二氧化碳吸收量，则二氧化碳浓度为组内因子。均衡设计：各个组观测数目相同组内因子：时间是两水平（五周、六周）的组内因子，单因素组内方差分析：每个患者在所有水平下都进行了测量。重复测量方差分析：每个受试者不止一次被测量。多元方差分析：因变量不止一个多元协方差分析：协变量存在。协变量：白氏抑郁症量表（BDI）记录了他们的抑郁水平，那么你可以评测疗法类型的影响前，对任何抑郁水平的组间差异进行统计性调整。协方差分析：BDI称为协变量，该设计为协方差分析。ANOVA（）模型拟合ANOVA和回归方法都是独立发展而来，但是从函数形式上看，他们都是广义线性模型的特例。lm()函数也能分析ANOVA模。但是一般用aov()函数展示结果。aov(formula，data=dataframe)表达式中各项的顺序因子不止一个，并且是非平衡设计；存在协变量。出现任意一种情况，等式右边的变量都及其他变量相关，此时，无法清晰的划分他们对因变量的影响。表达式的顺序很重要#

y

~

A

#

y

~

x

+

A

#

y

~

A

*

B

#

y

~

x1

+

x2

+

A*B

#

y

~

B

+

A

#

y

~

A

+

Error(Subject/A)

#

y

~

B

*

W

+

Error(Subject/W)对双因素方差分析，若观测数不同，则Y~A*B及Y~B*A的结果不同。1序贯型Y~A+B+A：BR中的ANOVA表的结果将评价：A对Y的影响控制A时，B对Y的影响。控制A和B的主效应时，A及B的交互效应。2分层型：A根据B调整，B依据A调整，A:B交互项同时根据A和B调整3边界型：A根据B和A:B做调整，A:B交互项根据A和B调整样本大小越不平衡，效应项的顺序对结果的影响越大。一般来说越基础的效应越需要放在表达式前面。具体而言，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互项，以此类推。对主效应，越基础性的变量越应该放在表达式前面，因此性别放在处理方式之前。单因素方差分析比较分类因子定义的两个或者多个组别中因变量的均值。例子：50个患者接受5种治疗方案，求哪种药物疗法降低胆固醇（响应变量）最多绘制带有置信区间的均值图形。绘制带有95%置信区间的各疗法均值。>library(gplots)>plotmeans(response~trt,xlab="Treatment",ylab="Response",main="MeanPlot\nwith95%CI")>detach(cholesterol)95%置信区间的各疗法均值药物E降低胆固醇最多，而1time降低最少多重比较F检验表明5种药物疗效不同，但是并没有说明哪种疗效及其他疗效不同。TukeyHSD(fit)

par(las=2)

par(mar=c(5,8,4,2))

plot(TukeyHSD(fit))疑问：用什么函数得出组间差异置信区间的，Ｆ检验？方法二：#有相同字母的组说明均值差异不显著library(multcomp)par(mar=c(5,4,6,2))tuk<-glht(fit,linfct=mcp(trt="Tukey"))plot(cld(tuk,level=.05),col="lightgrey")评估检验的假设条件单因素方差分析结果准确性，依赖于可做单因素方差分析的假设。假设有假设因变量服从正态分布各组方差相等具体检验方法假设因变量服从正态分布用QQ图检验正态性各组方差相等用方差齐次性检验函数（对离群点敏感）用方差齐次性检验函数有Bartlett检验、Fligner-Killeen检验、Brown-Forsythe检验（HH包中hov()函数）。选一个就行，结果相同。要看看是否有离群点使得各组间方差没有显著不同故真的是各组方差数据没有显著差异。结论：该数据可以用ANOVA模型拟合得非常好。单因素协方差分析包含了一个或者多个定量的协变量。举个例子就是：降雨量（t）=K*温度（t）+e其中，t是自变量时间，降雨量（t）是因变量，而温度（t）则是协变量K为一个常数。例子：自变量（dose）怀孕小组分为四组，每组接受（0,5,50或者500）的药物处理。产下幼仔重量均值为因变量，怀孕时间为协变量。例子：自变量（dose）怀孕小组分为四组，每组接受（0,5,50或者500）的药物处理。产下幼仔重量均值为因变量，怀孕时间为协变量。通过F检验，检验两组的方差是否有显著性差异。S^2=∑(X-X平均)^2/(n-1)两组数据就能得到两个S^2值，S大^2和S小^2F=S大^2/S小^2由表中f大和f小（f为自由度n-1),查得F表，然后计算的F值及查表得到的F表值比较，如果F<F表表明两组数据没有显著差异；F≥F表表明两组数据存在显著差异置信度95%时F值(单边)f大f小2345678910∞2345678910∞19.09.556.945.795.144.744.464.264.103.0019.169.286.595.414.764.354.073.863.713.6019.259.126.395.194.534.123.843.633.482.3719.309.016.265.054.393.973.693.483.333.2119.338.946.164.954.283.873.583.373.222.1019.368.886.094.884.213.793.503.293.142.0119.378.846.044.824.513.733.443.233.071.9419.388.816.004.784.103.683.393.183.021.8819.398.785.964.744.063.633.343.132.971.8319.58.535.634.363.673.232.932.712.541.00获取调整的各组均值，即去除协变量效应后各组均值。F检验表明了不同剂量下各幼鼠体重均值不同，但是没有告诉我们哪种处理方式及其他方式不同。（此处是下多大剂量及其他组不同）T值为正值（2.581）说明未用药比用药体重高。其他对照调节rbind()详见help(glht)GeneralLinearHypotheses一般性线性假设平估检验的假设条件假设正态性同方差性回归斜率相同：本例中假定四个处理组通过怀孕时间（协变量）来预测出生体重（因变量）的回归斜率相同。自变量一般为类别型变量。（协变量定义时就有这一点，协变量及因变量必须呈现线性关系）ANCOVA模型包含怀孕时间*剂量的交互项时如：（Weight=gesttime+dose+gesttime:dose），可对回归斜率的同质性进行检验。交互效应若显著，则意味着怀孕时间和幼崽出生体重间的关系依赖于药物剂量的水平。当然不显著才行···交互不显著，则模型约为：Weight=gesttime+dose此时回归斜率就基本相同，只是截距不同。若假设不成立：可以尝试变换协变量或因变量或使用能对每个斜率独立解释的模型或使用不需要假设回归斜率同质性的非参数ancova方法。结果可视化检验斜率同质性的可视化方法install.packages("HH")library(HH)ancova(weight~gesttime+dose,data=litter)双因素方差分析例子：60只豚鼠，分别采用两种喂食方法（橙汁及维C），各种喂食方法中抗坏血酸含量有三种水平（0.5mg\1mg\2mg）,每种处理方式分配给10只豚鼠。因变量为牙齿长度。主效应和交互效应都非常显著意味着什么？即喂食方式、喂食剂量、以及二者交互都影响方差一致性。是吗？可视化处理：带置信区间的均值以上只显示了交互效应，未显示主效应。推荐下图方法：library(HH)interaction2wt(len~supp*dose)模型的假设检验及均值比较及上节类似。重复测量方差分析受试者不止被测量一次。一般情况含一个组间因子，一个组内因子。par(las=2)par(mar=c(10,4,4,2))with(w1b1,interaction.plot(conc,Type,uptake,type="b",col=c("red","blue"),pch=c(16,18),main="InteractionPlotforPlantTypeandConcentration"))boxplot(uptake~Type*conc,data=w1b1,col=c("gold","green"),main="ChilledQuebecandMississippiPlants",ylab="Carbondioxideuptakerate(umol/m^2sec)")二氧化碳浓度及植物类型对二氧化碳吸收的交互效应解释：ChilledQuebec比MississippiPlants二氧化碳吸收能力强，随着浓度增加，这种差异越来越明显。多元方差分析当因变量不止一个时，可用多元方差分析。单因素多元方差分析例子：研究鼓舞中卡路里、脂肪、糖是否会因为储存位置（1代表底层货架，2中层，3高层）不同而发生变化。卡路里、脂肪、糖为因变量，储存位置为自变量。上面仅说了不同货架上含卡路里脂肪糖均是不同的。可以补充个均值比较（如TukeyHSD）来判断对于每个因变量，哪个货架及其他货架是不同的。评估假设检验单因素多元方差分析假设：多元正态性向量x及均值u.x及均值的马氏距离平方服从自由度为p的卡方分布该例子中。P=3QQ图展示卡方分布的分位数，横纵坐标分别是样本量及马氏距离平方值。方差—协方差矩阵同质性：各组的协方差矩阵相同。通常可用Box’sM来评估该假设。但R中么有。对多元正态性敏感，很容易拒绝该假设。检验多元离群点library(mvoutlier)outliers<-aq.plot(y)outliers稳健多元方差分析如果多元正态性或者方差—协方差均值假设都不满足，又担心多元离群点，可以用稳健或者非参数版的MANOVA检验。不知道具体用啥模型用回归做ANOVA没看十、功效分析主要内容：判断所需样本量计算效应值评价统计功效样本大小：对于研究有X个受试者，值得做吗？到底需要多少个受试者显著性：H0为真，拒绝H0的概率（1型错误）功效：H0假设为假时，将H0拒绝的概率：（1-11型错误）效应值：依据所用假设检验不同，如：t检验，效应值为d=假设检验概述功效分析T检验方差分析相关性线性模型比例检验卡方检验效应值功效分析图形十一、中级绘图1散点图散点图矩阵高密度散点图（重合点很多）用封箱、颜色、透明度来指明任意点上重叠点的数目封箱图密度通过颜色标识三维散点图气泡图一般通过三维图展示三个定量变量之间关系，还有一种思路：先创建一个二维散点图，然后用点的大小来代表第三个变量的值，这便是气泡图。Plot()调用即添加新图Line()在原图形上添加信息绘图参数乱了，重启就恢复默认了2折线图#代码清单11-4展示五种橘树随时间推移的生长状况的折线图Orange$Tree<-as.numeric(Orange$Tree)ntrees<-max(Orange$Tree)xrange<-range(Orange$age)yrange<-range(Orange$circumference)plot(xrange,yrange,type="n",xlab="Age(days)",ylab="Circumference(mm)")colors<-rainbow(ntrees)linetype<-c(1:ntrees)plotchar<-seq(18,18+ntrees,1)for(iin1:ntrees){tree<-subset(Orange,Tree==i)lines(tree$age,tree$circumference,type="b",lwd=2,lty=linetype[i],col=colors[i],pch=plotchar[i])title("TreeGrowth","exampleoflineplot")legend(xrange[1],yrange[2],1:ntrees,cex=0.8,col=colors,pch=plotchar,lty=linetype,title="Tree"比较五种橘树生长状况的折现图3相关图顺