R语言方法总结.docx
《R语言方法总结.docx》由会员分享,可在线阅读,更多相关《R语言方法总结.docx(21页珍藏版)》请在冰豆网上搜索。
R语言方法总结
计算描述性统计量:
1、summary():
例:
summary(mtcars[vars])
summary()函数提供了最小值、最大值、四分位数和数值型变量的均值,以及因子向量和逻
辑型向量的频数统计。
2、apply()函数或sapply()函数
计算所选择的任意描述性统计量。
mean、sd、var、min、max、median、length、range
和quantile。
函数fivenum()可返回图基五数总括(Tukey’sfive-numbersummary,即最小值、
下四分位数、中位数、上四分位数和最大值)。
sapply()
例:
mystats<-function(x,na.omit=FALSE){
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
kurt<-sum((x-m)^4/s^4)/n-3
return(c(n=n,mean=m,stdev=s,skew=skew,kurtosis=kurt))
}
sapply(mtcars[vars],mystats)
3、describe():
Hmisc包:
返回变量和观测的数量、缺失值和唯一值的数目、平均值、
分位数,以及五个最大的值和五个最小的值。
例:
library(Hmisc)
describe(mtcars[vars])
4、stat.desc():
pastecs包
若basic=TRUE(默认值),则计算其中所有值、空值、缺失值的数量,以及最小值、最大值、值域,还有总和。
若desc=TRUE(同样也是默认值),则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数。
若norm=TRUE(不是默认的),则返回正态分布统计量,包括偏度和峰度(以及它们的统计显著程度)和Shapiro–Wilk正态检验结果。
这里使用了p值来计算平均数的置信区间(默认置信度为0.95:
例:
library(pastecs)
stat.desc(mtcars[vars])
5、describe():
psych包
计算非缺失值的数量、平均数、标准差、中位数、截尾均值、绝对中位差、最小值、最大值、值域、偏度、峰度和平均值的标准误
例:
library(psych)
describe(mtcars[vars])
分组计算描述性统计量
1、aggregate():
例:
aggregate(mtcars[vars],by=list(am=mtcars$am),mean)
2、by():
例:
dstats<-function(x)(c(mean=mean(x),sd=sd(x)))
by(mtcars[vars],mtcars$am,dstats)
by(mtcars[,vars],mtcars$am,plyr:
:
colwis(dstats))
3、summaryBy():
doBy包
例library(doBy)
summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)
4、describe.by():
doBy包(describe.by()函数不允许指定任意函数,)
例:
library(psych)
describe.by(mtcars[vars],mtcars$am)
5、reshape包分组:
(重铸和融合)
例:
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)
频数表和列联表
1、table():
生成简单的频数统计表
mytable<-with(Arthritis,table(Improved))
Mytable
2、prop.table():
频数转化为比例值
prop.table(mytable)
3、prop.table()*100:
转化为百分比
prop.table(mytable)*100
二维列联表
4、table(A,B)/xtabs(~A+b,data=mydata)
例:
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
5、margin.table()和prop.table():
函数分别生成边际频数和比例
(1:
行,2:
列)
行和与行比例
margin.table(mytable,1)
prop.table(mytable,1)
列和与列比例
margin.table(mytable,2)
prop.table(mytable,2)
prop.table(mytable)
6、addmargins():
函数为这些表格添加边际和
addmargins(mytable)
admargins(prop.table(mytable))
addmargins(prop.table(mytable,1),2)
addmargins(prop.table(mytable,2,1)
7.crossTable():
gmodels包
例:
library(gmodels)
CrossTable(Arthritis$Treatment,Arthritis$Improved)
多维列联表
1、table()和xtabs():
都可以基于三个或更多的类别型变量生成多维列联表。
2、ftable():
例:
mytable<-xtabs(~Treatment+Sex+Improved,data=Arthritis)
mytable
ftable(mytable)
margin.table(mytable,1)
margin.table(mytable,2)
margin.table(mytable,3)
margin.table(mytable,c(1,3))
ftable(prop.table(mytable,c(1,2)))
ftable(addmargins(prop.table(mytable,c(1,2)),3))
gtable(addmargins(prop.table(mytable,c(1,2)),3))*100
独立检验
1、卡方独立性检验:
chisq.test()
例:
library(vcd)
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
chisq.test(mytable)
mytable<-xtabs(~Improved+Sex,data=Arthritis)
chisq.test(mytable)
2、Fisher精确检验:
fisher.test()
例:
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
fisher.test(mytable)
3、Cochran-Mantel—Haenszel检验:
mantelhaen.test()
例:
mytable<-xtabs(~Treatment+Improved+Sex,data=Arthritis)
mantelhaen.test(mytable)
相关性度量
1、assocstats():
例:
library(vcd)
mytable<-xtabs(~Treatment+Improved,data=Arthritis)
assocstats(mytable)
2、cor():
函数可以计算这三种相关系数,
3、cov():
函数可用来计算协方差
例:
states<-state.x77[,1:
6]
cov(states)
cor(states)
cor(states,method="spearman")
x<-states[,c("Population","Income","Illiteracy","HSGrad")]
y<-states[,c("LifeExp","Murder")]
cor(x,y)
4、pcor():
偏相关ggm包
例:
library(ggm)
pcor(c(1,5,2,3,6),cov(states))
相关性的显著性检验
1、cor.test()
其中的x和y为要检验相关性的变量,alternative则用来指定进行双侧检验或单侧检验(取值为"two.side"、"less"或"greater"),而method用以指定要计算的相关类型("pearson"、
"kendall"或"spearman")当研究的假设为总体的相关系数小于0时,请使用alternative=
"less"。
在研究的假设为总体的相关系数大于0时,应使用alternative="greater"。
在默认情况下,假设为alternative="two.side"(总体相关系数不等于0)。
例:
cor.test(states[,3],states[,5])
2、corr.test():
可以为Pearson、Spearman或Kendall相关计算相关矩阵和显著性水平。
例:
library(psych)
corr.test(states,use="complete")
3、pcor.test():
psych包
t检验
1、t.test(y~x,data)(独立样本)
例:
library(MASS)
t.test(Prob~So,data=UScrime)
2、t.test(y1,y2,paired=TRUE)(非独立)
例:
library(MASS)
sapply(UScrime[c("U1","U2")],function(x)(c(mean=mean(x),
sd=sd(x))))
with(UScrime,t.test(U1,U2,paired=TRUE))
组间差异的非参数检验
两组的比较:
1、wilcox.test(y~x,data):
评估观测是否是从相同的概率分布中抽得
例:
with(UScrime,by(Prob,So,median))
wilcox.test(Prob~So,data=UScrime)
2、wilcox.test(y1,y2,paried=TRUE):
它适用于两组成对数据和无法保证正态性假设的情境。
例:
sapply(UScrime[c("U1","U2")],median)
with(UScrime,wilcox.test(U1,U2,paired=TRUE))
多于两组的比较:
1、kruskal.test(y~A,data):
各组独立
例:
states<-as.data.frame(cbind(state.region,state.x77))
kruskal.test(Illiteracy~state.region,data=states)
2、friedman.test(y~A|B,data):
各组不独立
非参数多组比较:
1、npmc():
npmc包
例:
class<-state.region
var<-state.x77[,c("Illiteracy")]
mydata<-as.data.frame(cbind(class,var))
rm(class,var)
library(npmc)
summary(npmc(mydata),type="BF")
aggregate(mydata,by=list(mydata$class),median)
回归
用一个或多个预测变量(也称自变量或解释变量)来预测响应变量(也称因变量、效标变量或结果变量)的方法。
1、lm():
拟合回归模型lm(y~x1+x2+x3,data)
简单线性回归
1、lm():
(data是数据框)
例:
fit<-lm(weight~height,data=women)
summary(fit)
women$weight
fitted(fit)
residuals(fit)
plot(women$height,women$weight,main="WomenAge30-39",
xlab="Height(ininches)",ylab="Weight(inpounds)")
多项式回归
例:
fit2<-lm(weight~height+I(height^2),data=women)
summary(fit2)
plot(women$height,women$weight,main="WomenAge30-39",
xlab="Height(ininches)",ylab="Weight(inlbs)")
lines(women$height,fitted(fit2))
2、scatterplot():
绘制二元关系图
例:
library(car)
scatterplot(weight~height,data=women,spread=FALSE,
lty.smooth=2,pch=19,main="WomenAge30-39",xlab="Height(inches)",
ylab="Weight(lbs.)")
多元线性回归
1、scatterplotMatrix():
car包
scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图,并添加平滑(loess)
和线性拟合曲线。
对角线区域绘制每个变量的密度图和轴须图。
例:
fit<-lm(Murder~Population+Illiteracy+Income+
Frost,data=states)
有交互项的多元线性回归
例:
fit<-lm(mpg~hp+wt+hp:
wt,data=mtcars)
summary(fit)
1、effect():
effects包:
展示交互项的结果
term即模型要画的项,mod为通过lm()拟合的模型,xlevels是一个列表,指定变量要设定的常量值,multiline=TRUE选项表示添加相应直线。
例:
library(effects)
plot(effect("hp:
wt",fit,xlevels=list(wt=c(2.2,3.2,4.2))),
multiline=TRUE)
回归诊断
1、confint():
求模型参数的置信区间
例:
fit<-lm(Murder~Population+Illiteracy+Income+
Frost,data=states)
confint(fit)
2、plot():
生成评价模型拟合情况的图形
例:
fit<-lm(weight~height,data=women)
par(mfrow=c(2,2))
plot(fit)
3、lm():
删除观测点
例:
newfit<-lm(weight~height+I(height^2),data=women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)
par(opar)
gvlma包提供了对所有线性模型假设进行检验的方法
检验正态性:
4、qqPlot():
car包:
学生化残差(studentizedresidual,也称学生化删除残差或折叠化残差)
例:
library(car)
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
qqPlot(fit,labels=row.names(states),id.method="identify",simulate=TRUE,main="Q-QPlot")
注:
id.method="identify"选项能够交互式绘图
5、fitted():
提取模型的拟合值
例:
fitted(fit)[“Nevada”]
6、residuals():
二项式回归模型的残差
例:
residuals(fit)[“Nevada”]
7、residplot():
生成学生化残差柱状图(即直方图),并添加正态曲线、核密度曲线和轴须图。
它不需要加载car包
例:
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)
误差的独立性
8、durbinWatsonTest():
验证独立性
例:
durbinWatsonTest(fit)
验证线性
9、crPlots():
car包成分残差图也称偏残差图
例:
crPlots(fit)
同方差性(car包的两个函数)
10、ncvTest():
生成一个计分检验,零假设为误差方差不变,备择假设为误差方差随着拟合值水平的变化而变化。
若检验显著,则说明存在异方差性
11、spreadLevelPlot():
添加了最佳拟合曲线的散点图,展示标准化残差绝对值与拟合值的关系。
例:
library(car)
ncvTest(fit)
spreadLevelPlot(fit)
线性模型假设的综合验证
1、gvlma():
gvlma包:
线性模型假设进行综合验证,同时还能做偏斜度、峰度和异方差性的评价
例:
library(gvlma)
gvmodel<-gvlma(fit)
summary(gvmodel)
多重共线性
1、vif():
car包:
函数提供VIF值,
>2就表明存在多重共线性问题
例:
vif(fit)
sqrt(vif(fit))>2
异常观测值
1、outlierTest():
car包:
求得最大标准化残差绝对值Bonferroni调整后的p值
例:
library(car)
outlierTest(fit)
高杠杆值点
1、hat.plot():
观测点的帽子值大于帽子均值的2或3倍,即可以认定为高杠杆值点
例:
hat.plot<-function(fit){
p<-length(coefficients(fit))
n<-length(fitted(fit))
plot(hatvalues(fit),main="IndexPlotofHatValues")
abline(h=c(2,3)*p/n,col="red",lty=2)
identify(1:
n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(fit)
强影响点:
Cook’sD值大于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")
1、influencePlot():
car包:
离群点、杠杆值和强影响点的信息整合到一幅图形中
例:
influencePlot(fit,id.method="identify",main="InfluencePlot",
sub="CirclesizeisproportialtoCook'sDistance")
纵坐标超过+2或小于-2的州可被认为是离群点,
水平轴超过0.2或0.3的州有高杠杆值(通常为预测值的组合)。
圆圈大小与影响成比例,圆圈很大的点可能是对模型参数的估计造成的不成比例影响的强影响点
变量变换
1、powerTransform():
car包:
函数通过λ的最大似然估计来正态化变量
。
例:
library(car)
summary(powerTransform(states$Murder))
2、boxTidwell():
car包:
通过获得预测变量幂数的最大似然估计来改善线性关系
例:
library(car)
boxTidwell(Murder~Population+Illiteracy,data=states)
模型比较
1、anova():
基础包:
比较两个嵌套模型的拟合优度
例:
fit1<-lm(Murder~Population+Illiteracy+Income+
Frost,data=states)
fit2<-lm(Murder~Population+Illiteracy,data=states)
anova(fit2,fit1)
2、AIC():
AIC值越小的模型(可以不嵌套)要优先选择,它说明模型用较少的参数获得了足够的拟合度。
例:
fit1<-lm(Murder~Population+Illiteracy+Income+
Frost,data=states)
fit2<-lm(Murder~Population