OLS回歸是通過預測變量的加權和來預測量化的因變量,，其中權重是通過數(shù)據(jù)估計而得以的參數(shù)。

使殘差平方和最小

為能夠恰當?shù)亟忉孫LS模型的系數(shù),，數(shù)據(jù)必須滿足以下統(tǒng)計假設：

（1）      正態(tài)性對于固定的自變量,，因變量值成正態(tài)分布

（2）      獨立性 Yi值之間相互獨立

（3）      線性   因變量與自變量之間為線性相關

（4）      同方差性因變量的方差不隨自變量的水平不同而變化，即不變方差或同方差性

3. 用lm（）擬合回歸模型

[plain] view plain copy

1. fit<-lm(weight~height,data=women)  
2. summary(fit)  

[plain] view plain copy

1. fitted(fit)#擬合模型的預測值  

[plain] view plain copy

1. residuals(fit)#擬合模型的殘差值  

[plain] view plain copy

1. plot(women$height,women$weight,  
2.      xlab="Height （in inches）",  
3.      ylab="Weight（in pounds）")  
4. abline(fit)  

[plain] view plain copy

1. fit2<-lm(weight~height+I(height^2),data=women)  
2. summary(fit2)  

[plain] view plain copy

1. plot(women$height,women$weight,  
2.      xlab="Height（in inches）",  
3.      ylab="Weight（in lbs）")  
4. lines(women$height,fitted(fit2))  

[plain] view plain copy

1. scatterplot(weight~height,  
2.             data=women,  
3.             spread=FALSE,  
4.             lty.smooth=2,  
5.             pch=19,  
6.             main="Women Age 30-39",  
7.             xlab="Height (inches)",  
8.             ylab="Weight(lbs.)")  

[plain] view plain copy

1. states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])  

[plain] view plain copy

1. cor(states)  

[plain] view plain copy

1. library(car)  
2. scatterplotMatrix(states,spread=FALSE,lty.smooth=2,main="Scatter Plot Matrix")  

[plain] view plain copy

1. fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
2. summary(fit)  

[plain] view plain copy

1. fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)  
2. summary(fit)  

[plain] view plain copy

1. fit<-lm(mpg~hp+wt+hp:wt,data=mtcars)  
2. summary(fit)  
3.   
4. install.packages("effects")  
5. library(effects)  
6. plot(effect("hp:wt",fit,  
7.             list(wt=c(2.2,3.2,4.2))),multiline=TRUE)  

[plain] view plain copy

1. 簡單線性回歸  

[plain] view plain copy

1. fit<-lm(weight~height,data=women)  
2. par(mfrow=c(2,2))  
3. plot(fit)  

[plain] view plain copy

1. fit2<-lm(weight~height+I(height^2),data=women)  
2. par(mfrow=c(2,2))  
3. plot(fit2)  

[plain] view plain copy

1. library(car)  
2. fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
3. qqPlot(fit,labels=row.names(states),id.method="identify",simulate=TRUE,main="Q-Q Plot")  

[plain] view plain copy

1. residplot<-function(fit,nbreaks=10){  
2.   z<-rstudent(fit)  
3.   hist(z,breaks=nbreaks,freq=FALSE,  
4.        xlab="Studnetized Residual",  
5.        main="Distribution of Errors")  
6.   rug(jitter(z),col="brown")  
7.   curve(dnorm(x,mean=mean(z),sd=sd(z)),  
8.         add=TRUE,col="blue",lwd=2)  
9.   lines(density(z)$x,density(z)$y,  
10.         col="red",lwd=2,lty=2)  
11.   legend("topright",  
12.          legend=c("Normal Curve","Kernel Density Curve"),  
13.          lty=1:2,col=c("blue","red"),cex=0.7)}  
14. residplot(fit)  

[plain] view plain copy

1. durbinWatsonTest(fit)  

[plain] view plain copy

1. library(car)  
2. crPlots(fit)  

[plain] view plain copy

1. library(car)  
2. ncvTest(fit)  
3. spreadLevelPlot(fit)  

[plain] view plain copy

1. install.packages("gvlma")  
2. library(gvlma)  
3. gvmodel<-gvlma(fit)  
4. summary(gvmodel)  

[plain] view plain copy

1. library(car)  
2. vif(fit)  
3. sqrt(vif(fit))>2  

[plain] view plain copy

1. library(car)  
2. outlierTest(fit)  

[plain] view plain copy

1. hat.plot<-function(fit){  
2.   p<-length(coefficients(fit))  
3.   n<-length(fitted(fit))  
4.   plot(hatvalues(fit),main="Index Plot of Hat Values")  
5.   abline(h=c(2,3)*p/n,col="red",lty=2)  
6.   identify(1:n,hatvalues(fit),names(hatvalues(fit)))  
7. }  
8. hat.plot(fit)  

[plain] view plain copy

1. cutoff<-4/(nrow(states)-length(fit$coefficients)-2)  
2. plot(fit,which=4,cook.levels=cutoff)  
3. abline(h=cutoff,lty=2,col="red")  

[plain] view plain copy

1. library(car)  
2. avPlots(fit,ask=FALSE,onepage=TRUE,id.method="identify")  

[plain] view plain copy

1. library(car)  
2. influencePlot(fit,id.method="identify",main="Influence Plot",  
3.               sub="Circle size if proportional to Cook's distance")  

[plain] view plain copy

1. library(car)  
2. summary(powerTransform(states$Murder))  

[plain] view plain copy

1. fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
2. fit2<-lm(Murder~Population+Illiteracy,data=states)  
3. anova(fit2,fit1)  

[plain] view plain copy

1. fit1<-lm(Murder~Population+Illiteracy+Income+Frost,  
2.          data=states)  
3. fit2<-lm(Murder~Population+Illiteracy,data=states)  
4. AIC(fit1,fit2)  

[plain] view plain copy

1. 后向回歸  
2. library(MASS)  
3. fit1<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
4. stepAIC(fit,direction="backward")  

[plain] view plain copy

1. install.packages("leaps")  
2. library(leaps)  
3. leaps<-regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)  
4. plot(leaps,scale="adjr2")  

[plain] view plain copy

1. library(car)  
2. subsets(leaps,statistic="cp",main="Cp Plot for All Subsets Regression")  
3. abline(1,1,lty=2,col="red")  

[plain] view plain copy

1. install.packages("bootstrap")  
2. library(bootstrap)  
3. shrinkage<-function(fit,k=10){  
4.   require(bootstrap)  
5.   theta.fit<-function(x,y){lsfit(x,y)}  
6.   theta.predict<-function(fit,x){cbind(1,x)%*%fit$coef}  
7.   x<-fit$model[,2:ncol(fit$model)]  
8.   y<-fit$model[,1]  
9.   results<-crossval(x,y,theta.fit,theta.predict,ngroup=k)  
10.   r2<-cor(y,fit$fitted.values)^2  
11.   r2cv<-cor(y,results$cv.fit)^2  
12.   cat("Original R-square=",r2,"\n")  
13.   cat(k,"Fold Cross-Validated R-square=",r2cv,"\n")  
14.   cat("Change=",r2-r2cv,"\n")  
15. }  
16. fit<-lm(Murder~Population+Income+Illiteracy+Frost,data=states)  
17. shrinkage(fit)  
18. fit2<-lm(Murder~Population+Illiteracy,data=states)  
19. shrinkage(fit2)  

[plain] view plain copy

1. zstates<-as.data.frame(scale(states))  
2. zfit<-lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)  
3. coef(zfit)  

[plain] view plain copy

1. relweights<-function(fit,……){  
2.   R<-cor(fit$model)  
3.   nvar<-ncol(R)  
4.   rxx<-R[2:nvar,2:nvar]  
5.   rxy<-R[2:nvar,1]  
6.   svd<-eigen(rxx)  
7.   evec<-svd$vectors  
8.   ev<-svd$values  
9.   delta<-diag(sqrt(ev))  
10.   lambda<-evec%*%delta%*%t(evec)  
11.   lambdasq<-lambda^2  
12.   beta<-solve(lambda)%*%rxy  
13.   rsqrare<-colSums(beta^2)  
14.   rawwgt<-lambdasq%*%beta^2  
15.   import<-(rawwgt/rsquare)*100  
16.   lbls<-names(fit$model[2:nvar])  
17.   rownames(import)<-lbls  
18.   colnames(import)<-"Weight"  
19.   barplot(t(import),names.arg=lbls,  
20.           ylab="% of R-Square",  
21.           xlab="Predictor Variables",  
22.           main="Relative Importance of Predictor Variables",  
23.           sub=paste("R-Square=",round(rsquare,digits=3)),……)  
24.   return(import)  
25. }  
26. fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)  
27. relweights(fit,col="lightgrey")