R语言-回归分析相关函数

lm(formula,data)#formula是模型公式，data为数据框
#y~1+x1+x2和y~x1+x2是一样的，若想拟合齐次的，改成y~0+x1+x2或y~-1+x1+x2
summary(object)#object为lm函数生成的对象


#例1.
blood<-data.frame(
x1=c(76,91.5,85.5,82.5,79,80.5,74.5,79,85,76.5,82,95,92.5),
x2=c(50,20,20,30,30,50,60,50,40,55,40,40,20),
y=c(120,141,124,126,117,125,123,125,132,123,132,155,147))
lm.sol<-lm(y~1+x1+x2,data=blood)
summary(lm.sol)
#Call:
#lm(formula = y ~ 1 + x1 + x2, data = blood)


#Residuals:
#    Min      1Q  Median      3Q     Max 
#-4.0404 -1.0183  0.4640  0.6908  4.3274 


#Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
#(Intercept) -62.96336   16.99976  -3.704 0.004083 ** 
#x1            2.13656    0.17534  12.185 2.53e-07 ***
#x2            0.40022    0.08321   4.810 0.000713 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


#Residual standard error: 2.854 on 10 degrees of freedom
#Multiple R-squared:  0.9461,    Adjusted R-squared:  0.9354 
#F-statistic: 87.84 on 2 and 10 DF,  p-value: 4.531e-07
#***表示极为显著，**表示高度显著，*表示显著，.表示不太显著，没有标记表示不显著
#回归系数和回归方程都是显著的，所以最后的回归方程为：Y=-62.96+2.136X1+0.4002X2


predict(object,newdata,interval=c("none",confidence","prediction"),level=0.95)
#object为lm函数生成的对象,newdata为数据框，默认值为已知数据


#例2.（接着例1）
newdata<-data.frame(x1=80,x2=40)
predict(lm.sol,newdata,interval="prediction")
predict(lm.sol,newdata,interval="confidence")


anova()#计算方差分析表
coef()#提取系数
deciance()#计算残差平方和
formula()#提取模型公式
labels()#标记
proj()#给出数据在线性模型上的投影


#回归诊断
#例3.
data(anscombe)
ff<-y~x
for(i in 1:4){
 ff[2:3]<-lapply(paste(c("y","x"),i,sep=""),as.name)
 assign(paste("lm.",i,sep=""),lmi<-lm(ff,data=anscombe))
}
getcoef<-function(obj) summary(get(obj))$coef 
lapply(objects(pat="lm\\.[1-4]$"),getcoef)
#可以看到，四组数据的计算结果基本上相同
#事实上这四组数据完全不同，所以全部用线性模型是不合适的


#例4.（例3的绘图程序）
opar<-par(mfrow=c(2,2),mar=0.1+c(4,4,1,1),oma=c(0,0,2,0))
for(i in 1:4){
 ff[2:3]<-lapply(paste(c("y","x"),i,sep=""),as.name)
 plot(ff,data=anscombe,col="red",pch=21,bg="orange",cex=1.2,xlim=c(3,19),ylim=c(3,13))
 abline(get(paste("lm.",i,sep="")),col="blue")
}
mtext("Anscombe's 4 Regression data sets",outer=T,cex=1.5)
par(opar)


#残差检验：检查误差是否满足正态性和方差齐性
residuals(object)#简写resid(),计算回归模型残差
rstandard(model)#计算回归模型的标准化(也称内学生化)残差
rstudent(model)#计算回归模型的外学生化残差
#所谓外学生化残差就是删除第i个样本数据后得到的标准化残差
#object和model为lm生成的对象


#例5.（数据与例3相同）
attach(anscombe)
lm.sol1<-lm(y1~x1)
pre<-fitted.values(lm.sol1)
res<-residuals(lm.sol1)
rst<-rstandard(lm.sol1)
opar<-par(mfrow=c(1,2))
plot(pre,res,xlab="Fitted Values",ylab="Residuals")
plot(pre,rst,xlab="Fitted Values",ylab="Standardized Residuals")
par(opar)
detach(anscombe)


#影响分析
influence.measures(model)#model为lm生成的对象
dffits(model,infl,res)#DFFITS准则
dfbeta(model)
dfbetas(model)
covratio(model)#COVRATIO准则
cooks.distance(model)#cook距离
hatvalues(model)#帽子值
hat(x,intercept=T)#帽子矩阵，其中x为设计阵


#例6.（数据与例1相同）
blood<-data.frame(
x1=c(76,91.5,85.5,82.5,79,80.5,74.5,79,85,76.5,82,95,92.5),
x2=c(50,20,20,30,30,50,60,50,40,55,40,40,20),
y=c(120,141,124,126,117,125,123,125,132,123,132,155,147))
lm.sol<-lm(y~1+x1+x2,data=blood)
influence.measures(lm.sol)
opar<-par(mfrow=c(3,2),mar=0.4+c(4,4,1,1),oma=c(0,0,2,0))
plot(lm.sol,1:6)
par(opar)


#Box-Cox变换
boxcox(object,lambda=seq(-2,2,1/10),plotit=T)#object为lm函数生成的对象


#例7.
x<-c(679,292,1012,493,582,1156,997,2189,1097,2078,
1818,1700,747,2030,1643,414,354,1276,745,435,
540,874,1543,1029,710,1434,837,1748,1381,1428,
1255,1777,370,2316,1130,463,770,724,808,790,
783,406,1242,658,1746,468,1114,413,1787,3560,1495,2221,1526)
y<-c(0.79,0.44,0.56,0.79,2.7,3.64,4.73,9.5,5.34,6.85,
5.84,5.21,3.25,4.43,3.16,0.5,0.17,1.88,0.77,1.39,
0.56,1.56,5.28,0.64,4,0.31,4.2,4.88,3.48,7.58,
2.63,4.99,0.59,8.19,4.79,0.51,1.74,4.1,3.94,0.96,
3.29,0.44,3.24,2.14,5.71,0.64,1.9,0.51,8.33,14.94,5.11,3.85,3.93)
lm.sol2<-lm(y~x)
opar<-par(mfrow=c(2,2),mar=0.4+c(4,4,1,1),oma=c(0,0,2,0))
plot(fitted(lm.sol2),col="red")
library(MASS)
boxcox(lm.sol2,lambda=seq(0,1,by=0.1))
lambda1<-0.55 #从图中可以看出
y1<-(y^lambda1-1)/lambda1
lm.lam<-lm(y1~x)
plot(fitted(lm.lam),col="orange")
beta0<-lm.lam$coef[1]
beta1<-lm.lam$coef[2]
curve((1+lambda1*(beta0+beta1*x))^(1/lambda1),from=min(x),to=max(x),col="blue")
points(x,y,col="green")
mtext("Box-Cox",outer=T)
par(opar)


#多重共线性
#当设计矩阵存在多重共线性时，最小二乘估计的性质不够理想，有时甚至很坏
#这时需要新的估计方法，岭估计是最有影响且应用较为广泛的估计方法
lm.ridge(formula,data,lambda=0)#formula为模型公式，lambda为岭参数，默认为0(即最小二乘估计)


#例8.
france<-data.frame(
x1=c(149.3,161.2,171.5,175.5,180.8,190.7,202.1,212.4,226.1,231.9,239),
x2=c(4.2,4.1,3.1,3.1,1.1,2.2,2.1,5.6,5,5.1,0.7),
x3=c(108.1,114.8,123.2,126.9,132.1,137.7,146,154.1,162.3,164.3,167.6),
y=c(15.9,16.4,19,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3)
)
lm.rid<-lm.ridge(y~x1+x2+x3,data=france,lambda=seq(0,0.2,length=50))
par(mai=c(0.9,0.9,0.2,0.2))
plot(lm.rid)#从图中可以看出k=0.04时能满足岭迹选择k的条件
#确定参数的方法HKB、LW
select(lm.rid)#对于本问题，这两种方法得到的参数并不理想
lm.rid1<-lm.ridge(y~x1+x2+x3,data=france,lambda=0.04);lm.rid1