柚子快報(bào)邀請碼778899分享：R語言實(shí)戰(zhàn)-第八章回歸

http://yzkb.51969.com/

#第八章 回歸

#簡單線性回歸 #用到基礎(chǔ)包中的women數(shù)據(jù)集，研究身高與體重的關(guān)系 head(women) fit <- lm(weight~height,data=women) summary(fit) fitted(fit)#列出擬合模型的預(yù)測值 residuals(fit)#列出擬合模型的殘差值 plot(women$height,women$weight) abline(fit) #或者lines(women$height,fitted(fit)) 兩種方法的差異在右上角 #擬合后方差解釋率為99.1% print(fit)

#多項(xiàng)式回歸 觀察到women散點(diǎn)圖有曲線的趨勢，嘗試用多項(xiàng)式回歸看回歸結(jié)果怎樣 head(women) fit1 <- lm(weight~height+I(height^2),data=women) summary(fit1) fitted(fit1)#列出擬合模型的預(yù)測值 residuals(fit1)#列出擬合模型的殘差值 plot(women$height,women$weight) lines(women$height,fitted(fit1))#注意這里和簡單線性回歸繪制的方式不一樣 #擬合后方差解釋率為99.9%

#多元線性回歸 #以state.x77數(shù)據(jù)為例 研究人口 文盲率等對謀殺率的貢獻(xiàn)作用(未考慮預(yù)測變量的交互項(xiàng)) states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) #先檢查變量之間的相關(guān)性 cor(states) library(car) car::scatterplotMatrix(states,spread= FALSE,smoother.args = list(lty=2), ? ? ? ? ? ? ? ? ? ? ? ?main = "Scatterplot matrix") #非對角線繪制變量間的散點(diǎn)圖，并添加平滑和線性擬合曲線 #對角線區(qū)域繪制每個(gè)變量的密度圖和軸須圖

#使用lm()進(jìn)行多元線性回歸模型 fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) summary(fit) #forest/income不與murder呈線性關(guān)系，總體看來， #所有預(yù)測變量解釋了各州謀殺率57%的方差

#有交互項(xiàng)的多元線性回歸 #以mtcars數(shù)據(jù)集為例，探討馬力hp和車重wt對每英里耗油量mpg的影響 fit <- lm(mpg~hp+wt+hp:wt,data = mtcars) summary(fit) #交互項(xiàng)也很顯著，表明響應(yīng)變量與其中一個(gè)預(yù)測變量的關(guān)系依賴于另一個(gè)預(yù)測變量的水平 #effects::effect()展示交互項(xiàng)的結(jié)果 #plot(effect(term,mod,,xlevels),multiline=TRUE) library(effects) plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline=TRUE)

#回歸診斷 #基本的方法 confint states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) summary(fit) confint(fit) #標(biāo)準(zhǔn)方法 plot函數(shù) 檢驗(yàn)正態(tài)性 獨(dú)立性 線性和同方差性 fit <- lm(weight~height,data = women) par(mfrow=c(2,2)) plot(fit)

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

#刪除點(diǎn)后(刪除數(shù)據(jù)需謹(jǐn)慎，本應(yīng)該是模型去匹配數(shù)據(jù)，而不是反過來) newfit2 <- lm(weight~height+I(height^2),data = women[-c(13,15),]) par(mfrow=c(2,2)) plot(newfit2)

states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit2 <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) par(mfrow=c(2,2)) plot(fit2)

#改進(jìn)的方法? #1 正態(tài)性：car包qqPlot函數(shù) ##car包已經(jīng)更新原來的id.method="identify"已經(jīng)不可用了，變成了id=list() library(car) states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) #用下面的新代碼 car::qqPlot(fit,main="Q-Q Plot",id=list(method="identify", ? ? ? ?labels=row.names(states)),simulate=TRUE,pch=16)

#這里出現(xiàn)的“警告: 沒有點(diǎn)在0.25英尺內(nèi)”是Rstudio的問題， #放入R本來的環(huán)境運(yùn)行不存在該問題。建議進(jìn)入原環(huán)境運(yùn)行

#car::qqPlot(fit,labels=row.names(states),id.method="identify", # ? ? ? simulate=TRUE,main="Q-Q plot")無法交互作用

#Nevada在置信區(qū)間外，關(guān)注一下它： #觀察實(shí)際謀殺率和模擬的謀殺率差別 states["Nevada",] fitted(fit)["Nevada"] residuals(fit)["Nevada"] rstudent(fit)["Nevada"]

#利用residplot可視化誤差，生成學(xué)生化殘差柱狀圖 residplot <- function(fit,nbreaks=10){ ? z <- rstudent(fit) ? hist(z,breaks=nbreaks,freq=FALSE, ? ? ? xlab="studentized residual", ? ? ? main="distribution of errors") ? 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("normal curve","kernel density curve"), ? ? ? ? ?lty=1:2,col = c("blue","red"),cex=0.7) } residplot(fit)

#2 誤差的獨(dú)立性 最好的方法是依據(jù)數(shù)據(jù)收集方式的先驗(yàn)知識(shí)? #car::durbinWatsonTest也可以檢測相關(guān)性 library(car) car::durbinWatsonTest(fit) #如果加上simulate=TRUE則運(yùn)行結(jié)果與沒有時(shí)有些許不同 #p值不顯著 說明無自相關(guān)性

#3 線性通過成分殘差圖觀察 car::crPlots(fit) #若圖像存在非線性 則說明預(yù)測變量的函數(shù)形式可能建模不夠充分，可能需要添加一些曲線成分 #比如多項(xiàng)式 或?qū)σ粋€(gè)或多個(gè)變量進(jìn)行變換

#4 同方差性 car包中的ncvTest和spreadLevelPlot函數(shù) car::ncvTest(fit) #p值不顯著說明滿足同方差性 car::spreadLevelPlot(fit) #會(huì)顯示建議的變換，此處異方差性很不明顯因此建議冪次接近1，不需要進(jìn)行變換； #注意建議冪變換為0則使用對數(shù)變換

#線性模型假設(shè)的綜合驗(yàn)證 #gvlma::gvlma會(huì)給出綜合驗(yàn)證，并進(jìn)行偏斜度 峰度 異方差性的評價(jià) library(gvlma) gvmodel <- gvlma::gvlma(fit) summary(gvmodel)

#多重共線性? #car::vif 一般原則下根號(hào)下vif>2表明存在多重共線性問題 car::vif(fit) sqrt(car::vif(fit))>2 #均為FALSE 說明不存在多重共線性問題

#異常值觀測 #1 離群點(diǎn)? #前面用Q-Qplot判斷? #也可以用標(biāo)準(zhǔn)化殘差值絕對值>2認(rèn)為是離群點(diǎn)粗略判斷 #更好的方法 car::outlierTest car::outlierTest(fit) #Nevada的Bonferroni p被認(rèn)為是離群點(diǎn) 刪除后還需要在檢驗(yàn)是否有其他離群點(diǎn)存在

#高杠桿值點(diǎn) :由許多異常的預(yù)測變量值組合起來的 與響應(yīng)變量值沒有多大關(guān)系 #使用帽子統(tǒng)計(jì)量判斷 帽子均值為p/n 一般認(rèn)為觀測點(diǎn)的帽子值大于帽子均值的2或3倍就可以認(rèn)為是高杠桿值點(diǎn) hat.plot <- function(fit){ ? p <- length(coefficients(fit)) #模型估計(jì)的參數(shù)數(shù)目(包括截距) ? n <- length(fitted(fit)) #樣本量 ? plot(hatvalues(fit),main = "Index plot of hat values") ? abline(h=c(2,3)*p/n,col="red",lty=2) ? identify(1:n,hatvalues(fit),names(hatvalues(fit))) } hat.plot(fit) #高杠桿值點(diǎn)可能是強(qiáng)影響點(diǎn) 也可能不是 這要看他們是否是離群點(diǎn) #警告: 沒有點(diǎn)在0.25英尺內(nèi) 是Rstudio的問題，放入R本來的環(huán)境運(yùn)行不存在該問題 #點(diǎn)擊finish結(jié)束吧

#強(qiáng)影響點(diǎn) #Cook's D值大于4/(n-k-1)表示有強(qiáng)影響點(diǎn) #可以通過變量添加圖判斷

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)#length(fit$coefficients)包括了截距 所以這里-2 plot(fit,which = 4,cook,levels=cutoff) abline(h=cutoff,lty=2,col="red")

car::avPlots(fit,ask=FALSE,id.method="identify")

#利用car::influencePlot將離群點(diǎn) 杠桿值和強(qiáng)影響點(diǎn)信息整合到一幅圖 car::influencePlot(fit,id.method="identify") #縱坐標(biāo)在±2之外可被認(rèn)為是離群點(diǎn)；水平軸超過0.2或0.3具有高杠桿值； #圓圈大小與影響成比例，圓圈很大的點(diǎn)可能是模型參數(shù)的估計(jì)造成的不成比例影響的強(qiáng)影響點(diǎn)

#違背回歸假設(shè)的數(shù)據(jù)能改進(jìn)的措施 #刪除觀測點(diǎn) 變量變換 添加或刪除變量 使用其他回歸方法 #刪除觀測點(diǎn) ?持謹(jǐn)慎態(tài)度

#變量變換? library(car) states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) #用下面的新代碼 car::qqPlot(fit,main="Q-Q Plot",id=list(method="identify", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? labels=row.names(states)),simulate=TRUE,pch=16) #1 當(dāng)模型違反正態(tài)假設(shè)時(shí) 通?？梢詫憫?yīng)變量嘗試某種變換 car::powerTransform #2 當(dāng)模型違反線性假設(shè)時(shí) 通常可以對預(yù)測變量嘗試某種變換 car::boxTidwell

#1 響應(yīng)變量變換 Box-Cox正態(tài)變換 car::powerTransform library(car) summary(powerTransform(states$Murder)) #結(jié)果建議用y^0.6進(jìn)行變換 λ=0.6接近0.5 可以嘗試平方根轉(zhuǎn)換 #但LR test, lambda = (1) ?p=0.14512 即λ=1時(shí)也無法拒絕 證明本例可以不進(jìn)行變換

#2 預(yù)測變量變換 car::boxTidwell #用人口和文盲率預(yù)測謀殺率 car::boxTidwell(Murder~Population+Illiteracy,data=states) #結(jié)果建議population^0.87和illiteracy^1.36可以改善線性關(guān)系 #但計(jì)分檢驗(yàn)p值表明不需要進(jìn)行變換 #變量變換應(yīng)該謹(jǐn)慎 因?yàn)樽儞Q有意義才好解釋

#增刪變量 #刪除某個(gè)存在多重共線性的變量（某個(gè)變量平方根vif >2） #或者使用嶺回歸——多元回歸的變體 專門用來處理多重共線性問題

#嘗試其它回歸方法 #存在離群點(diǎn)或強(qiáng)影響點(diǎn):OLS回歸→穩(wěn)健回歸模型 #違背正態(tài)性假設(shè):使用非參數(shù)回歸模型 #存在顯著的非線性:嘗試非線性回歸模型 #違背誤差獨(dú)立性假設(shè):利用專門研究誤差結(jié)構(gòu)的模型 如時(shí)間序列模型或多層次回歸模型 #轉(zhuǎn)向廣泛應(yīng)用的廣義線性模型 它適用于許多OLS回歸假設(shè)不成立的情況 #多重共線性問題 使用嶺回歸——多元回歸的變體?

#選擇"最佳"的回歸模型 #模型比較 基礎(chǔ)包anova() AIC()函數(shù) #anova()需要嵌套模型的模型才可以作比較 states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) fit2 <- lm(Murder~Population+Illiteracy,data = states) anova(fit2,fit1) #檢驗(yàn)不顯著 表明可以刪除Income+Frost

#AIC()函數(shù) 不需要嵌套模型也可以比較 AIC(fit2,fit1) #AIC值較小的模型要優(yōu)先選擇

#變量選擇 #逐步回歸模型(向前 向后 向前向后) 全自己回歸 #逐步回歸模型 MASS::stepAIC() library(MASS) states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) MASS::stepAIC(fit,direction="backward") #逐步回歸存在爭議 因?yàn)椴皇敲恳粋€(gè)可能的模型都被評價(jià)了

#全子集回歸 #leaps::regsubsets()展示結(jié)果 library(leaps) leapss <- leaps::regsubsets(Murder~Population+Illiteracy+Income+Frost, ? ? ? ? ? ? ? ? ? ? ? ? ? ? data = states,nbest=4) plot(leapss,scale="adjr2") #雙預(yù)測變量模型(最上面的0.55)population和illiteracy是最佳模型

#car::subsets()展示結(jié)果 car::subsets(leapss,statistic="cp", ? ? ? ? ? ? ?main="Cp plot for all subsets regression") abline(1,1,tyl=2,col="red") #越好的模型里截距和斜率均為1的直線越近

#注 大部分情況中 全子集回歸都優(yōu)于逐步回歸，但是全子集回歸較慢 #我們需要認(rèn)識(shí)到擬合效果家而沒有意義的模型對我們毫無幫助，我們應(yīng)該借助自己對知識(shí)背景的理解才能獲得最理想的模型

#深層次分析? #交叉驗(yàn)證 評價(jià)回歸方程的泛化能力 bootstrap::crossval library(bootstrap) #生成shrinkage函數(shù)用以評價(jià)泛化能力 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("Chang =",r2-r2cv,"\n") } states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states) shrinkage(fit) fit2 <- lm(Murder~Population+Illiteracy,data = states) shrinkage(fit2) #可見 fit2的模型具有更好的泛化能力，Chang值更小(因?yàn)橛^測值被隨機(jī)分配到k個(gè)群組， #所以每次運(yùn)shrinkage函數(shù)結(jié)果有少許不同)?

#相對重要性 可以比較標(biāo)準(zhǔn)化回歸系數(shù)或相對權(quán)重 #比較標(biāo)準(zhǔn)化回歸系數(shù) scale(注scale函數(shù)返回的是矩陣 而lm函數(shù)要求數(shù)據(jù)框 需要進(jìn)行轉(zhuǎn)換) states <- as.data.frame(state.x77[,c("Murder","Population", ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?"Illiteracy","Income","Frost")]) zstates <- as.data.frame(scale(states)) zfit <- lm(Murder~Population+Illiteracy+Income+Frost,data = zstates) coef(zfit) # Illiteracy 標(biāo)準(zhǔn)化回歸系數(shù)最大 因此解釋率最高

#相對權(quán)重 #生成relweights函數(shù)用以預(yù)測變量的相對權(quán)重 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 ? import <- as.data.frame(import) ? row.names(import) <- names(fit$model[2:nvar]) ? names(import) <- "Weights" ? import <- import[order(import),1,drop=FALSE] ? dotchart(import$Weights,labels=row.names(import), ? ? ? ? ? xlab="% of R-Square",pch=19,#決定系數(shù) ? ? ? ? ? main="Relative importance of predictor variables",#預(yù)測變量的相對重要性 ? ? ? ? ? sub=paste("Total R-Square=",round(rsquare,digits = 3)),#總決定系數(shù) ? ? ? ? ? ...) ? return(import) } relweights(fit,col="blue") ?

柚子快報(bào)邀請碼778899分享：R語言實(shí)戰(zhàn)-第八章回歸

http://yzkb.51969.com/

文章鏈接