R实战笔记..

R实战笔记（转载）
1.R的使用
1.3.基础命令
函数说明
example(‘foo’)查看函数示例
data() 查看已加载可用的数据集
vignette(‘foo’)查看主题为foo可用的vignette文档getwd() 显示当前工作目录
setwd() 修改工作目录
ls() 显示工作空间的对象
rm(list=ls()) 删除全部的对象
options() 显示或设置当前选项
history(n) 显示最近的n个命令
save.image(‘myfile’)保存工作空间到myfile中，默认是.RData load(‘myfile’)读取一个工作空间到当前会话
source(‘file_name’)在当前会话中执行一个脚本
sink(‘file_name’)重定向输出到文件，append是否追加,split=T 则同时输出到屏幕和文件中
object.size(x)/1000000 查看变量占用的内存空间，M srirage.mode(x) 改变变量的存储类型
storage.mode(x）<- “integer”改为整数型，可以看到该对象的大小会变为原来的一半
memory.size() 查看现在的work space的内存使用
memory.limit() 查看系统规定的内存使用上限,注意，在32位的R中，封顶上限为4G，无法在一个程序上使用超过4G （数位上限）。

这种时候，可以考虑使用64位的版本
rm(object) 删除变量
gc() 做Garbage collection，否则内存是不会自动释放的，相当于你没做rm.
rm(list=ls()) 删除全部变量
图像输出
pdf(‘filename.pdf’) #重定向到图像输出
png(‘filename.png’)
jpeg(‘filename.jpeg’)
dev.off() #将输出返回到终端
2.数据结构
2.2.2.矩阵
matrix(vector, nrow=n, ncol=m, byrow=TRUE, dianames=list(rname, cname))
2.2.4.数据框
数据框绑定
attach()
detach()
with()
attach(mtcars)
plot(mpg, wt)
mpgdetach(mtcars)
with(mtcars, {
summary(mpg, disp, wt)
plot(mpg, disp)
plot(mpg, wt)
})
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10 2.2.6.列表的引用
list[[n]]
4.数据管理
1.data.frame(col1,col2…,stringAsFactors=FALSE)
创建数据框的时候，不要将字符串转成因子
2.创建新变量
mydata = transform(mydata,
sumx = x1+x2,
meanx = (x1+x2)/2)
library(dplyr)
mutate(mydata, sumx = x1+x2, meanx = (x1+x2)/2)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
4.3.变量重编码
!x：非x
x|y：或
x&y：与，和
isTRUE(x)：测试x是否为true
mydata$age[mydata$age == 99] <-NA
mydata$agecat[mydata$age >=55 & mydata$age <=75] <- 'middle aged' mydata = with(mydata, {
agecat <- NA
agecat[age>75] <-'elder'
agecat[age>=55 & age <=75] <- 'middle aged'
agecat[age<55] <- 'young' })
#epicalc包的recode函数#看看有哪些因子水平，或者说不重复的取值
list = unique(data3$arpu)
#构建old 列表
old <-
c("(50,80)","(10,20]","0","(20,50]","(0,5]","(100,150]","(150,200]","(200, 300]","(5,10]","[80,100]","300以上" )
#构建new 列表
new <- c("低值","低值","低值","低值","低值","高值","高值","高值","低值","高值","高值")
library(epicalc)
recode(vars = "arpu",old,new,dataframe)
#recode 函数其实是用new替换原来的old，所以建议先复制一列，然后在复制的列上进行recode 操作
#R的cut函数
aaa <- c(1,2,3,4,5,2,3,4,5,6,7)
cut(aaa, 3, b = 4, ordered = TRUE)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
∙20
∙21
∙22
∙23
∙24
∙25
∙26
∙27
∙28
4.4.变量重命名
#1
fit(mydata) #交互式修改`
#2library(reshape)
mydata = rename(mydata,
c(old_name1 = 'new_name1', old_name2 = 'new_name2',...))
#3
names(mydata)[2]='new_name'
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
4.5.缺失值
删除含有缺失值的行
na.omit(mydata)
4.6.日期
符号含义示例
%d 01-31天
%a 缩写的星期Mon
%A 非缩写的星期Monday
%m 月份00-12
%b 缩写的月份Jan
%B 非缩写的月份January
%y 两位数的年份07
%Y 4位数的年份2007
as.Data(c('2007-06-22','2004-02-13'))
#系统日期
Sys.Date()
date()
format(Sys.Date(), format='%Y%m%d')
#日期加减
end_date - start_date
diff(today, that_day, units='weeks')
#lubridate包
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13 4.8.排序
mydata = mydata[order[mydata$age),]
attach(mydata)
newdata = mydata[order(gender, -age),]detach(mydata) ∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
4.9.合并数据集
# merge,一种内连接，innder join
total = merge(data1, data2, by=c('id','country'))
rbind()
cbind()
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
4.10.子集
mydata[, c(6:10)]
mydata[, c(var1,var2,...)]
#丢弃变量
vars = names(mydata) %in% c('q3','q4')
mydata[!vars]
mydata[c(-8,-9)]
#删除变量
mydata$q4 <- NULL
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
4.10.进入观测
mydata = mydata[,1:3]
mydata[ which(mydata$gender == 'M' & mydata$age>30),] attach(mydata)
mydata[which(gender == 'M' & age>30),]detach(mydata) ∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
subset(mydata,age>=35 | age<24, select=c(q1,q2,q3,q4)) subset(mydata, gender =='M' & age >25, select =gender:q4) ∙ 1
∙ 2
∙ 3
4.11.随机抽样
# 无放回随机抽样，抽10个样本
mydata[sample(1:nrow(mydata),10, replace=FALSE), ] ∙ 1
∙ 2
4.12.sql语句
library(sqldf)
newdata = sqldf('select * from mydata where carb=1 order by mpg', s=TRUE)
new_data = sqldf('select avg(mpg), avg(disp) as avg_disp, gear from mtcars where cyl in (4,6) group by gear')
∙ 1
∙ 2
∙ 3
∙ 4
5.高级数据管理
5.2.数值与字符处理函数
5.2.1.数学函数
函数描述
abs(x) 绝对值
sqrt(x) 平方根
x^2 平方
ceiling(x) 向上取整
floor(x) 向下取整
trunc(x) 向零取整
round(x, digits=n) 直接截断保留n位小数
signif(x) 四舍五入保留n位小数
cos(x),sin(x),tan(x) 余弦，正弦，正切
acos(x),asin(x),atan(x) 反余弦，反正弦，反正切
cosh(x),sinh(x),tanh(x) 双曲余弦，双曲正弦，双曲正切acosh(x),asinh(x),atanh(x) 反双曲余弦，反双曲正弦，反双曲正切
log(x,base=n) n为底的对数,log(x)为自然对数，log10(x)为10常用对数
exp(x) 指数
5.2.2.统计函数
函数描述mean(x) 平均数
median(x) 中位数
sd(x) 标准差
var(x) 方差
mad(x) 绝对中位差
quantile(x,probs=c(0.3,0.84)) 分位数
range(x) 范围
sum(x) 求和
diff(x,lag=n) 滞后差分
min(x),max(x) 最小最大值
scale(x,center=T, scale=T) 中心化(center=T)或标准化(center=T, scale=T)
其他均值和方差的标准化
scale(mydata)*sd+mean
5.2.3.概率函数
d=密度函数
p=分布函数
q=分位数函数
r=生产随机数(随机偏差)
分布名称缩写分布名称缩写正态分布norm 非中心卡方分布chisq Beta分布beta logistic分布logis
二项分布binom 多项分布multinom 柯西分布cauchy 负二项分布nbinom 指数分布exp 泊松分布pois
F分布 f Wilcoxon符号秩分布signrank Gamma分布gamma t分布t
几何分布geom 均匀分布unif
超几何分布hyper Weibull分布weibull
分布名称缩写分布名称缩写
对数正态分布lnorm Wilcoxon秩和分布wilcox
设定随机数种子
set.seed(1234)
5.2.4.字符串处理函
函数描述
nchar(x)
x的字符数量，注意和
length(x)的区别
substr(x,start,stop)
提取或替换子字符串，
substr(‘abcde’,2,4)<-‘222
22’
grep(pattern,x,ignore.case=F,fixed=F)
正则表达式搜索，返回值
为匹配的下标
sub(pattern,replacement,x,ignore.case=F,fixe
d=F)
搜索并替换,如果
fixed=T,则pattern是一个
文本字符串
strsplit(x,split,fixed=F)
在split处分割字符串，
返回一个列表，用unlist()
变成向量
paste(x,sep=”)连接字符串，指定分隔符toupper(x) 转大写
tolower(x) 转小写
5.2.5.其他实用函数
函数说明
length(x) 长度
seq(from,to,by) 生产序列
rep(x,n,each=n) 将x重复n次,指定each将会排序
cut(x,n) 分割成n水平的因子
pretty(x,n)创建美观的分割点，在绘图中常用
cat(…,file=’myfile’,append=T)|连接…中的对象，并将其输出到屏幕或者文件中
5.3.apply-by系列
5.3.1.apply
# 1=对每行操作，2=对每列操作
apply(mydata, margin=1, fun=myfunc)
∙ 1
∙ 2
5.3.2.by,分组汇总操作
# 针对1:4列，根据species因子分组，求均值
by(iris[,1:4],Species,mean)
∙ 1
∙ 2
pply
l=list(a=1:10,b=11:20)
lapply(l,mean)
∙ 1
∙ 2
5.3.4.sapply
lapply返回的是一个含有两个元素a和b的list，而sapply返回的是一个含有元素**“a”++和**“b”++的vector，或者列名为a和b的矩阵。

l=list(a=1:10,b=11:20)
l.mean=sapply(l,mean)
class(l.mean)#提取元素a的均值
l.mean[['a']]
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
5.3.5.tapply
attach(iris)#根据sprcies进行分类，计算petal的均值
tapply(iris$Petal.Length,Species,mean)detach(iris)
∙ 1
∙ 2
∙ 3
∙ 4
7基本统计分析
7.1.1描述性统计分析
# 1.查看一般的统计量
vars = c('mpg','hp','wt')
summary(mtcars[vars])
# 2.对一列或多列应用多个函数
myfunc <- function(x, na.omit=Flase){
if(na.omit)
x = x[!is.na(x)]
m=mean(x)
n=length(x)
s=sd(x)
skew=sum((x-m)^3/s^3)/n
return(c(n=n,mean=m,stdev=s,skew=skew))
sapply(mtcars[vars],myfunc)
# 3.其他包里面的描述性统计量library(Hmisc)
describe(mtcars[vars])
library(psych)
describe(mtcars[vars])
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
∙20
∙21
∙22
∙23
∙24
7.1.2分组汇总
vars = c('mpg','hp','wt')
# 1.简单的aggregate
aggregate(mtcars[vars], by=(am=mtcars$am,vs=mtcats$vs),mean) # 2.一般化的by
myfunc = function(x)(c(mean=mean(x),sd=sd(x)))
by(mtcars[vars], mtcars$am, myfunc)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
7.2.1列连表分析
library(vcd) mydata = Arthritis
# 1.一维列联表# 1.1 生成频数统计表
mytable = with(mydata, table(Improved))
# 1.2 将频数统计转成百分比
prop.table(mytable)
# 2.二维列联表# 2.1 table的用法
mytable = table(A,B)
# 2.2 xtable的用法
mytable = xtabs(~ A+B, data=dataframe)
mytable = xtabs(~Treatment+Imporved, data=mydata)
# 2.3 添加边际和
addmargins(mytable)
# 2.4 转成百分比
prop.table(mytable)
# 3.三维列联表
mytable = xtabs(~Treatment+Imporved+Sex, data=mydata)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
∙20
∙21
∙22
∙23
∙24
∙25
∙26
∙27
7.2.2 检验类别型(分类型变量)独立性的方法。

1.卡方独立性检验
2.Fisher精确检验
3.Cochran-Mantel-Heanszel卡方检验
# 1.卡方检验library(vcd)
mydata = Arthritis
# 1.2先转列联表
mytable = xtabs(~Treatment+Imporved, data=mydata)
chisq.test(mytable)
#p<0.05,说明不独立。

# 2.Fisher精确检验的原假设是：边界固定的列联表中行和列是相互独立的。

fisher.test(mytable)
#p<0.05,说明不独立。

# 3.Cochran-Mantel-Heanszel卡方检验的原假设是：两个名义变量在第三个变量的每一层中都是条件独立的。

mytable = xtabs(~Treatment+Imporved+Sex, data=mydata) mantelhean.test(mytable)
#p<0.05,说明不独立。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
∙20
7.2.3 相关性的度量
既然变量不独立，那么如何衡量相关性的强弱呢?
library(vcd)
mydata = Arthritis
# 先转列联表
mytable = xtabs(~Treatment+Imporved, data=mydata)
assocstats(mytable)
# phi系数，列联系数，Cramer's V系数，越大说明相关性较强。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
7.3 相关
描述定量变量之间的关系。

包括pearson系数，Spearman系数，Kendall 相关系数，偏相关系数，多分格(polychoric)系数和多系列（polyserial)相关系数。

7.3.1 pearson,spearman,kendall 相关
∙pearson：描述两个变量之间的线性相关程度
∙spearman：描述分级定序变量之间的相关程度
∙Kendall’s Tau：也是一种非参数的等级相关度量
mydata = state.x77[,1:6]
# 1.pearson相关
cor(mydata)
# 2.spearman相关
cor(mydata,method='spearman')
# 3.kendall相关
cor(mydata,method='kendall')
# 注意:得到相关系数后，还需要对相关系数进行独立性检验。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12 7.3.2 偏相关
指控制一个或多个变量，看另外两个变量之间的相互关系。

library(ggm)
pcor(c(1,5,2,3,6), cov(mydata))
#c(1,5,2,3,6)：前面两个是要计算的偏相关的变量的下标，剩余的是要控制(排除影响)的变量下标。

∙ 1
∙ 2
∙ 3
∙ 4
7.3.3相关性的显著性检验
常用的原假设是：变量间不相关（即总体的相关系数为0）。

cor.test(mydata[,3],mydata[,5])
#p<0.05，说明真的是存在相关性。

#如果总体相关系数小于0，加参数alternative='less'，大于0则alternative='greater'，默认是双侧检验(two.side)#参数method='pearson(默认),spearman,kendall'
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
# 一次计算全部变量的相关系数和显著性检验library(psych)
corr.test(mydata, use='complete')
# 得到correlation matrix(相关系数矩阵)，probility value(显著性P值矩阵)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
# 偏相关的检验library(psych)
corr_value = pcor(c(1,5,2,3,6), cov(mydata))
pcor.test(r,q,n)# r即corr_value# q是要控制的变量数(以数值表示位置)# n是样本的大小
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
7.4 t检验
研究中最常见的就是对两个组进行比较，包括用没有某种方法，药物，同一种方法前后的变化，两种不同药物的比较，两种不同工艺良品率的比较等。

这里关注结果变量为连续型的组间比较，并假设其呈现正态分布。

7.4.1 独立样本的T检验
针对两组独立样本T检验可以用于检验两个总体的均值是否相等。

原假设是均值相等，即这两种方法没有什么差别，两种药物的效果是一样的，前后没有明显的变化。

t.test(y~x,data)
y:数值型变量，x:二分变量
默认方差不等，可以用参数var.equal=T 假定方差相等。

library(MASS)
t.test(Prob~So, data=UScrime)
#p<0.05，说明南方和北方各州的监禁率是不等的。

∙ 1
∙ 2
∙ 3
∙ 4
7.4.2 非独立样本的T检验
比如年轻的比年长的失业率是否更高？这个就是不独立的。

在两组的观测之间相关时，获得的是一个非独立组设计，前-后测试，重复性测量设计同样会产生非独立的组。

假定组间的差异呈现正态分布。

t.test(y1,y2, paired=T)
y1,y2是两个非独立组的数值向量
p<0.05,说明二者的均值是不一样的，即有差别。

7.4.3 多余两组的比较
用方差分析(ANOVA).
7.5 组间差异的非参数检验
如果无法满足t检验或者方差分析的参数假设，就需要非参数方法。

7.5.1 两组的比较
∙ 1.若两组数据独立，用Wilcoxon秩和检验(即Mann-Whitney U检验)来评估是否从相同的概率分布中抽的，即检验是否
来自相同的总体。

# y是数值型变量，x是二分变量
wilcox.test(y~x, data)
#y1,y2是各组的结果变量
wilcox.test(y1,y2)
# p<0.05，说明来自不同的总体，即均值什么的是不一样的，二者有差异。

library(MASS)
wilcox.test(Prob~So, data=UScrime)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙ 2.wilcoxon符号秩和检验是非独立样本T检验的一种非参数替代方法，适用于两组成对数据和无法保证正态分布假
设的情景。

library(MASS)
with(UScrime, wilcox.test(U1,U2,paried=T))
∙ 1
∙ 2
7.5.2 多于两组的比较
比如要比较美国4个地区的文盲率。

这称为单向设计。

如果无法满足ANOVA设计的假设，那么可以使用非参数方法来评估组间的差异，如果各组独立，则Kruskal-Wallis检验，如果不独立(如重复测量设计或随机分组设计)，则Friedman检验。

Kruskal.test(y~A, data)
y是数值型结果变量，A是拥有两个及以上的分组变量
friedman.test(y~A|B, data)
y是数值型结果变量，A是分组变量，B是一个用以认定匹配观测的区组变量(blocking variable)
mydata = as.data.frame(cbind(state.region, state.x77))
kruskal.test(Illiteracy~state.region, data=mydata)
# p<0.5,说明四个地区的文盲率各不相同。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
但是没有告诉你哪些地区显著与其他地区不同。

这时候需要使用Mann-Whitney U检验每次比较两组数据。

一个更优雅的方法是：在控制了犯第一类错误的概率的前提下，执行可以同步进行的多组比较，这样可以直接完成所有组之间的成对比较。

使用npmc包。

class = state.region
var = state.x77[,c('Illiteracy')]
mydata = as.data.frame(cbind(class,var))
library(npmc)
summary(npmc(mydata),type='BF')
# p.value.2s双侧概率查看两两比较的结果，p<0.05说明存在明显差异。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
8回归
8.1回归的多面性
回归类
型
用途
简单线
性
用一个量化解释变量预测一个量化的响应变量
多项式
用一个量化解释变量预测一个量化的响应变量，模型的关系
是N阶多项式
多元线
性
用两个多多个量化的解释变量预测一个量化的响应变量
多变量用一个或多个解释变量预测多个响应变量
Logistic 用一个或多个解释变量预测一个类别型响应变量
泊松用一个或多个解释变量预测一个代表频数的响应变量
cox比例
风险
用一个或多个解释变量预测一个时间(死亡，失败或旧病复
发)发生的时间
时间序
列
对误差项相关的时间序列建模
非线性
用一个或多个量化解释变量预测一个量化的响应变量，模型
的关系非线性的
非参数
用一个或多个量化解释变量预测一个量化的响应变量，模型
源自数据形式，不能实现设定
稳健
用一个或多个量化解释变量预测一个量化的响应变量，能抵
御强影响点的干扰
8.2 最小二乘法(OLS)回归
Y^=β0^+β1^X1i+...++βi^Xki
求解使得残差平方和最小
∑(Yi−Y^i)2=∑(Yiβ0^+β1^X1i+...++βi^Xki)2=∑ε2
OLS的基本假设
- 正态性：对于固定的自变量值，因变量成正态分布。

- 独立性：Yi值(或残差)之间相互独立
- 线性：因变量和自变量之间为线性关系
- 同方差性：因变量的方差不随自变量的水平不同而变化，也称为不变方差。

对你和模型非常游泳的函数
函数说明
lm 拟合函数myfit=lm(formula, mydata)
summary(myfit)
展示模型的详细结果。

Coefficients为模型参数，
Multiple R-squared为拟合优度(R2),p值为模型整体
显著性检验值,P<0.05说明模型整体通过检验，参数
的检验在Coefficients中
coefficients(myfit) 模型的拟合参数
confint(myfit) 模型参数的置信区间
fitted(myfit) 模型的预测值
residuals(myfit) 模型的残差值
anova(myfit)
生成拟合模型的方差分析表，或者比较两个或者更
多拟合模型的方差分析表
vcov(myfit) 列出模型参数的协方差矩阵
AIC(myfit) 输出赤池统计量
plot(myfit) 生成拟合模型的诊断图
predict() 用模型去预测新的数据集
8.2.2简单的线性回归
fit = lm(weight~height, data=women)
summary(fit)
# 查看预测值
fitted(fit)
# 查看模型的残差
residuals(fit)
# 生成拟合模型的诊断图
plot(women$height, women$weight,xlab='x轴'，ylab='y轴')
# 添加拟合直线
abline(fit)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
8.2.3多项式回归
fit2 = lm(weight~height+I(height^2), data=women)
summary(fit2)
# 生成拟合模型的诊断图
plot(women$height, women$weight,xlab='x轴'，ylab='y轴')
# 添加拟合直线
lines(women$height,fitted(fit2))
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
8.2.4多元线性回归
最好先检查下变量之间的相关性
mydata =
as.data.frame(state.x77[,c('Murder','Population','Illiteracy','Income','Fros t')])
# 检查变量之间的相关性
cor(mydata)
# 变量相关关系的散点图矩阵library(car)
scatterplotMatrix(mydata,spread=F,lty.smooth=2,main='scatter plot matrix')
# 多元线性回归
fit = lm(Murder~Population+Illiteracy+Income+Frost, data=mydata) summary(fit)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
8.2.5有交互项的多元线性回归
fit = lm(mpg~hp+wt+hp:wt, data=mtcars) summary(fit)
∙ 1
∙ 2
∙ 3
回归诊断
8.3.1标准方法
最简单的就是使用plot(fit)绘制模型的拟合诊断图。

fit = lm(weigth~height, data=women)
par(mfrow=c(2,2))
plot(fit)
∙ 1
∙ 2
∙ 3
生成的四幅图对应的OLS假设
1. Residuals vs Fitted <–> 线性
2. Normal Q-Q <–> 正态性
3. Scale-Location <–> 同方差性
4. Residuals vs Leverage (残差与杠杆图，关注离群点)
1.正态性：当预测变量值固定时，因变量成正态分布，则残
差值也应该是一个均值为0的正态分布。

Q-Q图是在正态
分布对应的值下，标准化残差的概率图，Q-Q图上的点应
该落在45°直线上。

2.独立性：从上面的图是无法看出因变量是否相互独立，只
能从收集的数据集中验证。

3.线性：若y-x线性相关，那么残差值和预测值(拟合值)就没
有任何的系统关联。

除了白噪声，模型应该包含数据中的
所有方差。

Residuals vs Fitted图却有一个曲线关系，说明
应该对模型增加一个二次项。

4.同方差性：若满足同方差，则Scale-Location图中，水平线
周围的点应该随机分布。

8.3.2改进的方法
car包中的回归诊断函数
函数说明
qqPlot() 分位数比较图
durbinWatsonTest() 对误差自相关做Durbin-Watson检验crPots() 成分与残差图
ncvTest() 对非恒定的误差方法做得分分析spreadLevelPlot() 分散水平检验
outlierTest() Bonferroni离群点检测
avPlots() 添加的变量图像
inluencePlot() 回归影像图
scatterplot() 增强的散点图
scatterplotMatrix() 增强的散点图矩阵
vif() 方差膨胀因子
1.正态性
qqplot()画出n-p-1个自由度的t分布下的学生化残差，也成为学生化删除残差或折叠化残差。

n是样本量，p是回归参数的数目(包括截距项)。

library(car)
mydata =
as.data.frame(state.x77[,c('Murder','Population','Illiteracy','Income','Fros t')])
fit = lm(Murder~Population+Illiteracy+Income+Frost, data=mydata)
# 绘制95%置信区间的学生化残差图
qqplot(fit, labels=s(mydata),id.method='identity', simulate=T) # 在95%区间外的点，是模型高估或者低估该点(对应的现实事件是这个地点的谋杀率)。

mydata['Nevada']
# 查看模型预测值
fitted(fit)['Nevada']
# 查看残差
residuals(fit)['Nevada']
rstudent(fit)['Nevada']
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
2.误差的独立性
判断因变量值(或残差)是否相互独立，最好的方法是依据收集数据方式的先验知识。

Durbin-Watson检验，能检测误差的序列相关性。

durbinWatsonTest(fit)
# p<0.05,说明误差不独立。

lag=1表明数据集里面每个数据都是与其后一个数据进行比较的。

该检验适用于时间独立的数据，对于非聚集型的数据并不适用。

∙ 1
∙ 2
∙ 3
∙ 4
3.线性
通过成分残差图(也称偏残差图)，可以看到因变量和自变量之间是否呈现非线性关系。

也可以看看是否有不同于已设定线性模型的系统偏差。

library(car)
crPlots(fit)
#若图像存在非线性，说明可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分。

∙ 1
∙ 2
∙ 3
∙ 4
4.同方差性
ncvTest()函数生成一个计分检验，零假设是误差方差不变，备选假设是误差方差随着拟合值水平变化而变化。

spreadLevelPlot()创建一个添加了最佳拟合曲线的散点图，而且给出调整建议。

library(car)
ncvTest(fit)# p>0.05,说明方差相等，满足模型假设spreadLevelPlot(fit)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
8.3.3线性模型假设的综合验证
给出模型假设一个单独的综合检验(通过或者不通过)，如果不通过，需要哟个前面的方法去判断哪些假设没有被满足。

library(gvlma)
gvmodel=gvlma(fit)
summary(gvmodel)
# Global Stat的P值>0.05,说明满足OLS的所有假设。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
8.3.4多重共线性
使用方差膨胀因子vif,如果vif>4则说明存在多重共线性。

library(var)
vif(fit)
∙ 1
∙ 2
8.4.1离群值
1.方法1；在8.3.
2.1的qq图中，落在95%置信区间外的是
离群值。

2.方法2：标准化残差>2或<-2的也可能是离群值
library(car)
outlierTest(fit)
# Bonferonni P <0.05,说明是离群值。

∙ 1
∙ 2
∙ 3
∙ 4
8.4.2 高杠杆值点
即与其他预测变量有关的离群点，换句话说，他们是由许多异常的预测变量值组合起来，与响应变量值没有关系。

高杠杆值点是通过帽子统计量判断。

帽子均值是p/n,p是模型预计的参数数目(包括截矩)，n是样本量。

hat.plot = function(fit){
p=length(coefficients(fit))
n=length(fitted(fit))
plot(hatvalues(fit),main='index plot of hat values'
abline(h=c(2,3)*p/n, col='red',lty=2)
identity(1:n,hatvalues(fit),names(hatvalues(fit)))
}
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
8.4.3强影响点
即对模型参数估计值影响有些比例失衡的点，例如移除模型的一个观测点时模型会发生巨大的变化，那么就要检测下有没有强影响点的存在了。

1. 通过Cook距离，或称D统计量。

Cook’s D>4/(n-k-1)就是强影响点
2. 变量添加图
cutoff = 4/(nrow(mydata)-length(fiti$coefficients)-2)
plot(fit,which=4,cook.levels=cutoff)
abline(h=cutoff,lty=2,col='red')
∙ 1
∙ 2
∙ 3
将离群点，杠杆值，强影响点整合在一张图中。

library(car)
influencePlot(fit,id.method='identify')
# 纵坐标在[-2,2]之外的可认为是离群点，水平轴>0.2或0.3的有搞杠杆值。

圆圈的大小和影响成正比，圆圈很大的点可能是对模型参数的估计造成不成比例影响的强影响点。

∙ 1
∙ 2
∙ 3
∙ 4
8.6选择最佳的回归模型
8.6.1模型比较
AIC准则，越小越好，说明选用较少的参数获得足够好的拟合度。

fit1 = lm(Murder~Population+Illiteracy+Income+Frost, data=mydata)
fit2 = lm(Murder~Population+Illiteracy, data=mydata)
AIC(fit1,fit2)
∙ 1
∙ 2
∙ 3
∙ 4
8.6.2变量选择
1.逐步回归
library(MASS)
fit = lm(Murder~Population+Illiteracy+Income+Frost, data=mydata)
# 1.逐步回归
stepAIC(fit, direction='backward')
#然后根据给出的AIC变化钻则最佳的参数模型
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
2.全子集回归
library(leaps)
leaps=regsubsets(Murder~Population+Illiteracy+Income+Frost,
data=mydata, nbest=4)
# 根据调整R2挑选模型
plot(leaps,scale='adjust R2')
# 根据Mallows Cp统计量挑选模型# 越靠近y=x+1直线的越好。

但是从实际业务角度出发，有时候又不是这样的。

library(car)
subsets(leaps,statistic='cp')
abline(1,1,lty=2,col='red')
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
8.7交叉验证
交叉验证的函数
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]
resuts=corssval(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') #初始R方
cat(k,'fold cross-validated R-square=',r2cv,'\n') #交叉验证的R方}
# 1.先创建回归方差
fit = lm(Murder~Population+Illiteracy+Income+Frost, data=mydata)
# 2.交叉验证
shrinkage(fit,10)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
∙20
∙21
∙22
∙23
8.8相对重要性
1.标准化的回归系数
这是比较简单的方法。

在分析之前，先用scale()将数据标准化为mean=0,std=1的数据集，这样用R回归得到的就是标准化的回归系数。

（注意，scale得到的是矩阵，而lm接收的是数据库，所以要转换一下）
zmydata = as.data.frame(scale(states))
zfit = lm(Murder~Population+Illiteracy+Income+Frost, data=mydata) coef(fit)
∙ 1
∙ 2
∙ 3
2.相对权重法
它是对所有可能子模型添加一个预测变量引起的R2平均增加量的一个近似值。

# 1.自定义函数
relweights = functionc(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
delta=diag(sqrt(ev))
lambda=evec %*% delta %*% t(evec)
lambdasq=lambda^2
beta=solve(lambda) %*% rxy
rsquare=colSums(beta^2)
rawwgt=lambadasq %*% beta^2
import=(rawwgt/rsquare)*100
lbls=names(fit$model[2:nvar])
rownames(import)=lbls
rownames="weights"
barplot(t(import),names.arg=lbls,
ylab='% of R-square',
xlab='predctor variables'
main='relative importance of predictor variables',
sub=paste('R-square=',round(rsquare,digits=3)),
...)
return(import)
}
# 1.先拟合模型
fit=lm(Murder~Population+Illiteracy+Income+Frost, data=mydata)
# 2.返回变量的权重
relweights(fit,col='lightgrey')
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙17
∙18
∙19
∙20
∙21
∙22
∙23
∙24
∙25
∙26
∙27
∙28
∙29
∙30
∙31
∙32
9方差分析
当包含的因子是解释变量时，关注的重点僵尸组别差异的分析，这种分析法叫方差分析(ANOVA).
患者时间
5周
疗法CBT s1,s2,s3,s4,s5
EMD
R
s6,s7,s8,s9,s10
治疗方案是两水平(CBT,EMDR)的组间因子，之所以称作组间因子，是因为每个患者都仅分配到其中的一个方案中。

s是受试者(患者)。

因为每种治疗方案的观测数目都一样，所以也成为“均衡设计”，若观测数不同，则称为“非均衡设计”。

如果仅有一个类别型变量，则称为单因素方差分析。

疗法是组间因子，所以称为组间方差分析。

时间是两水平的组内因子，所以称为组内方差分析。

因为患者在不同的时间被测量，也称为重复测量方差分析。

疗法和时间组成的双因素方差分析。

若因子设计包括组内和组间因子，又称作混合模型方差分析。

疗法和时间都作为因子时，既可以分析疗法的影响(时间跨度平均)和时间的影响(疗法类型跨度的平均)，又可以分析疗法和时间的交互影响。

前两者叫主效应；后者叫交互效应。

方差分析主要通过F检验进行效果评测，原假设是组内或组间没有差异。

混淆因素：抑郁程度是接受治疗后的因变量，同时抑郁程度也会影响治疗的效果，所以抑郁程度也可以是组间差异，这时抑郁程度就是混淆因素。

干扰变数：如果你对混淆因素不感兴趣，那它就是干扰变数。

协变量：如果在评测疗法的影响前，对任何抑郁水平的组间差异进行统计调整。

这就是协变量。

协方差分析：协变量对应的方差分析。

如果因变量不止一个，就是多元方差分析(MANOVA),如果协变量也存在，就是多元协方差分析(MANCOVA)。

9.2ANOVA模型
aov(formula, data=dataframe)
设计formula
设计formula
单因素ANOVA y~A
含单个协变量的单因素ANCOVA y~x+A
双因素ANOVA y~A*B
含两个协变量的双因素ANCOVA y~x1+x2+A*B
随机化区组y~B+A(B是区组因子) 单因素组内ANOVA y~A+Error(Subject/A)
含单个组内因子(W)和单个组间因子(B)的重
复测量ANOVA
y~B*W+Error(Subject/B) 顺序很重要
y~A+B+A:B
1.序贯型：效应根据表达式中出现的效应作调整。

A不做调
整，B根据A调整，A:B交互项根据A和B调整。

2.分层型：效应根据同水平或低水平的效应做调整。

A根据
B调整，B根据A做调整，A:B同时根据A和B做调整。

3.边界型：每个效应根据其他各效应做相应的调整。

A根据
B和A:B做调整，A:B根据A和B做调整。

R默认采用序贯型方法，spss，sas默认使用边界型方法。

R中的ANOVA表的结果将评价：
1.A对y的影响
2.控制A时，B对y的影响
3.控制A和B时，A和B的交互效应
样本大小越不平衡，效应项的顺序对效果的影响越大。

一般来说，越基础的效应项越需要放在表达式的前面。

具体来讲，首先是协变量，然后是主效应，接着是双因素的交互项，再接着是三因素的交互效应，以此类推。

对于主效应，越基础的变量越应该放在表达式的前面，因此性别要放在处理方式之前。

有一个基本的准则：若研究设计不是正交的(也就是说，因子和或协变量相关)，一定要谨慎设置效应的顺序。

car包的Anova()函数(注意不是anova()函数)提供分层型和边界型的方法，而aov()使用序贯型方法。

如果要保持和SPSS,SAS的方法保持一致，可以使用Anova()函数。

9.3单因素方差分析
library(multcomp)attach(cholesterol)
#各样本的大小
table(trt)
#各组的均值
aggregate(response,by=list(trt),FUN=mean)
#各组的标准差
aggregate(response,by=list(trt),FUN=sd)
#检验组间差异
fit=aov(response~trt)
summary(fit)
#p<0.05,说明各组之间有差异
# 绘制各组均值及其置信区间的图像library(gplots)
plotmeans(response~trt, xlab='treatment',ylab='response',main='mean plot with 95% CI')
detach(cholesterol)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13
∙14
∙15
∙16
∙19
∙20
∙21
∙22
∙23
9.3.1 多重比较
aov()函数只说明各组之间有差异，不能说明哪种方法和其他方法有差别。

这时候要对各组均值差异的成对检验
TukeyHSD(fit)
∙ 1
9.3.2评估检验的假设条件
单因素方差分析的假设条件：
1. 假设因变量服从正态分布
2. 各组方差相等
正态分布的检验
library(car)
qqplot(lm(response~trt,data=cholesterol),simulate=T,main='Q-Q
plot',label=F)
#注意，qqplot()函数需要用lm()拟合。

#数据落在95%的置信区间内，说明满足正态分布假设。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
方差齐性检验
bartlett.test(response~trt, data=cholesterol)
#p>0.05,说明组间方差并没有显著不同。

∙ 1
∙ 2
∙ 3
方差齐性检验对离群点非常敏感，需要用前面的方法检测离群点。

library(car)
outlierTest(fit)
∙ 1
∙ 2
9.4 单因素协方差分析
data(litter,packag='multcomp')attach(litter)
table(dose)
aggregate(weight,by=list(dose),FUN=mean)
#注意：要将协变量放在表达式的前面
fit=aov(weight~gesttime+dose)
summary(fit)
#发现：a,怀孕时间和幼崽的出生体重相关，b,控制怀孕时间，药物剂量与出生体重有关。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
由于使用了协变量，可能想知道获取调整的组均值，即去除协变量(gesttime)效应后的组均值。

library(effects)
effect('dose',fit) 9.4.1组间均值的成对比较
library(multcomp)
contrast = rbind('no drug vs drug'=c(3,-1,-1,-1))
summary(glht(fit,linfct=mcp(dose=contrast)))
# p<0.05,说明未用药比其他用药条件下出生体重高的结论。

#c(3,-1,-1,-1)设定第一组和其他三组的均值比较，其他对组可用rbind()函数添加，详见help(glht)
∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
9.4.2评估检验的假设条件
同样要求正态性和同方差性假设外，还假定回归斜率相同
正态性和方差齐性可以用9.3的方法。

下面检测回归斜率的同质性。

ANCOVA模型包含怀孕时间*剂量的交互项时，可以回归斜率的同质性进行检验，交互项显著时，则意味着怀孕时间和出生体重间的关系依赖药物剂量的水平。

library(multcomp)
fit2=aov(weight~gesttime*dose, data=litter)
summary(fit2)
#看gesttime:dose项的p值，p>0.05，说明交互不显著，支持斜率相等的假设。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
若假设不成立，可以尝试变化协变量或因变量，或使用能对每个斜率独立解释的模型，或使用不需要假设回归斜率同质性的非参数ANCOVA方法。

9.5双因素方差分析
背景：60只老鼠，分几组，用2钟喂食方法，每种喂食方法又有3种不同的水平，看牙齿的长度。

library(ToothGrowth)
#均衡设计，因为每组的样本量是一样的
table(supp,dose)
#查看各组的均值和方差
aggregate(len,by=list(supp,dose),FUN=mean)
aggregate(len,by=list(supp,dose),FUN=sd)
#双因素方差分析
fit=aov(len~supp*dose)
#查看主效应和交互效应都非常显著p<0.05
summary(fit)
#可视化结果library(gplots)
plotmeans(len~interaction(supp,dose,sep=' '),
connect=list(c(1,3,5),c(2,4,6)),
col=c('red','darkgreen'),
main='interaction plot with 95% CIs',
xlab='treatment and dose combination')#图形展示了均值，误差棒(95%的置信区间)和样本量大小。

∙ 1
∙ 2
∙ 3
∙ 4
∙ 5
∙ 6
∙7
∙8
∙9
∙10
∙11
∙12
∙13。