一般線性回歸---完整過程

##線性回歸包括幾個方面：數(shù)據(jù)觀察，初步擬合蜕青，交互作用苟蹈，多重共線性，回歸診斷右核，擬合中出現(xiàn)的離群點(diǎn)慧脱，杠桿值，強(qiáng)影響點(diǎn)贺喝，刪除后菱鸥，重新擬合宗兼，模型比較--之后可能要預(yù)測，因此需要數(shù)據(jù)交叉驗(yàn)證分成2部分氮采，一部分?jǐn)M合殷绍，一部分驗(yàn)證。

數(shù)據(jù)為state.x77數(shù)據(jù)集---本次利用ols 最小二乘法=最小平方法直線回歸擬合分析

###線性回歸

##線性回歸包括幾個方面：數(shù)據(jù)觀察鹊漠，初步擬合主到，交互作用，多重共線性躯概，回歸診斷登钥，擬合中出現(xiàn)的離群點(diǎn)，杠桿值娶靡，強(qiáng)影響點(diǎn)牧牢，刪除后，重新擬合姿锭，模型比較--

##之后可能要預(yù)測塔鳍，因此需要數(shù)據(jù)交叉驗(yàn)證分成2部分，一部分?jǐn)M合艾凯，一部分驗(yàn)證献幔。

rm(list = ls())

gc()

library(car)

###本次利用ols 最小二乘法=最小平方法直線回歸擬合分析---數(shù)據(jù)來自于state.x77數(shù)據(jù)集

state.x77

data<-state.x77

data<-as.data.frame(data)? ###lm一定要data.frame格式

###先看看數(shù)據(jù)特點(diǎn)? car包的scatterplotMatrix()

data1<-data[,c("Murder","Population","Income","Illiteracy","Frost","Area")]

scatterplotMatrix(data1)

###線性擬合

fit<-lm(formula = Murder ~ Population +Income+Illiteracy+Frost +Area, data = data)

###多重共線性：sqrt(vif(fit))>2那一列 刪去

sqrt(vif(fit))

summary(fit)? ###Population? Illiteracy 具有與murder的線性關(guān)系

###看看之間是不是具有交互項(xiàng)

summary(lm(Murder~Population+Illiteracy+Population:Illiteracy,data=data))? ###沒有交互項(xiàng)

summary(lm(Murder~Population+Illiteracy,data=data))? ##公式為：Murder=1.652e+00+2.242e-04*Population+4.081e+00*Illiteracy

###殘差-回歸診斷 線性 正態(tài)性 獨(dú)立性 方差齊性

op<-par(mfrow=c(2,2))

plot(lm(Murder~Population+Illiteracy,data=data))? ## 沒有強(qiáng)影響點(diǎn)(沒有大的cook值)，Nevada>2(y值)為離群點(diǎn)趾诗，刪去離群點(diǎn)

par(op)

###刪去Nevada 重新擬合

rownames(data)

data2<-data[rownames(data)!="Nevada",]

summary(lm(Murder~Population+Illiteracy,data=data2))

op<-par(mfrow=c(2,2))

plot(lm(Murder~Population+Illiteracy,data=data2))? ## 沒有強(qiáng)影響點(diǎn)(沒有大的cook值)蜡感，Nevada>2(y值)為離群點(diǎn)，刪去離群點(diǎn)

par(op)

###改進(jìn)的殘差分析

class(data2)

str(data2)

fit_1<-lm(Murder~Population+Illiteracy,data=data2)

qqPlot(fit_1,id.method="identify",labels=row.names(data2)) ###qqplot正態(tài)性檢驗(yàn)

durbinWatsonTest(fit_1)? ##獨(dú)立性 durbin-watson檢驗(yàn) p>0.05,殘差獨(dú)立

crPlots(fit_1)? ###線性

ncvTest(fit_1)? ###方差齊性 p>0.05

spreadLevelPlot(fit)? ###方差齊性作圖 水平隨機(jī)分布的直線

###異常觀察值檢測 -總和方法

influencePlot(fit_1,id.method="identify")

###最佳回歸模型確定--全子集回歸

library(leaps)

leaps<-regsubsets(Murder ~ Population +Income+Illiteracy+Frost +Area, data = data,nbest = 5)

plot(leaps,scale="adjr2")? ##最好是2個自變量x

library(car)

subsets(leaps,statistic = "cp")

abline(1,1,lty=2,col="red")

###以上分析只是為了在數(shù)據(jù)基礎(chǔ)上建模和解釋參加模擬的數(shù)據(jù)恃泪；如果要模擬的曲線用來預(yù)測將要出來的數(shù)據(jù)時(shí)怎么辦郑兴？-深度分析之交叉驗(yàn)證

? ##交叉驗(yàn)證--9成數(shù)據(jù)來模擬，1成數(shù)據(jù)來驗(yàn)證

##set.seed()保證操作的可重復(fù)性贝乎，別人也用1234情连，產(chǎn)生的隨機(jī)數(shù)就和你的一樣了

set.seed(1234)

##選取訓(xùn)練集--從nrow（df）中即699個數(shù)字中，無放回（如果放回览效，replace=T）抽取0.7*nrow(df)個數(shù)字

train<-sample(nrow(data),0.7*nrow(data))

##提取出訓(xùn)練集的列表

df.train<-data[train,]

###提取驗(yàn)證集的列表

df.validate<-data[-train,]

##看看各個列表的數(shù)目

table(df.train$Murder)

table(df.validate$Murder)

##擬合：

fit_logit<-lm(Murder~.,data = df.train,family = binomial())

summary(fit_logit)

##用模擬的數(shù)據(jù)在新的數(shù)據(jù)集中却舀，進(jìn)行數(shù)據(jù)的驗(yàn)證: type = "response" 直接返回預(yù)測的概率值0~1之間

prob<-predict(fit_logit,df.validate,type="response")

prob

#3將概率大于0.5的定義為惡性腫瘤；概率小于0.5的定義為良性腫瘤,

logit.pred<-factor(prob>.5,levels = c(FALSE,TRUE),labels = c("benign","malignant"))

logit.pred

##3得出實(shí)際與預(yù)測的交叉表

logit.perf<-table(df.validate$class,logit.pred,dnn=c("Actual","Predicted"))

###預(yù)測出118個良锤灿，76個惡性

####? 準(zhǔn)確率為(76+118)/(129+81)=0.92 (76+118)/200=0.97

(76+118)/(129+81) ###有問題

(76+118)/200

###再回歸來看 有幾個模擬概率>0.05挽拔，不滿足，可以刪除再模擬但校，也可以用下面的方法

logit.fit.reduced<-step(fit_logit)

###哪個x值貢獻(xiàn)最大呢螃诅？scale先規(guī)范化為均值為1，方差為0的數(shù)據(jù)集

data_3<-scale(data)? ? ? ? ? ? #scale數(shù)據(jù)

data_3<-as.data.frame(data_3)

coef(lm(formula = Murder ~ Population +Income+Illiteracy+Frost +Area, data = data_3))

數(shù)據(jù)如下：

A B

25 7

25 9

25 9

27 12

27 14

27 16

24 16

30 14

30 16

31 16

30 17

31 19

30 21

28 24

32 15

32 16

32 17

32 25

34 27

34 15

34 15

35 15

35 16

34 19

35 18

36 17

37 18

38 20

40 22

39 25

43 24

?setwd("c:/Users/xx/Desktop/Homework_4A/")

> data<-read.table(file ="data_1.txt",header=T,sep = "\t")

> ##散點(diǎn)圖

> attach(data)

The following objects are masked from data (pos = 5):

? ? A, B

The following objects are masked from data (pos = 6):

? ? A, B

The following objects are masked from data (pos = 8):

? ? A, B

> ###pearson線性相關(guān)關(guān)系分析

> cor.test(data$A,data$B)

Pearson's product-moment correlation

data: ?data$A and data$B

t = 4.1427, df = 29, p-value = 0.0002712

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

?0.3257757 0.7927878

sample estimates:

? ? ? cor?

0.6097313?

> pic<-lm(B~A,data = data)

> print(pic)

Call:

lm(formula = B ~ A, data = data)

Coefficients:

(Intercept) ? ? ? ? ? ?A ?

? ? -2.3350 ? ? ? 0.6113 ?

> summary(pic)

Call:

lm(formula = B ~ A, data = data)

Residuals:

? ? Min ? ? ?1Q ?Median ? ? ?3Q ? ? Max?

-5.9469 -2.4765 -0.6145 ?1.4137 ?9.2193?

Coefficients:

? ? ? ? ? ? Estimate Std. Error t value Pr(>|t|) ? ?

(Intercept) ?-2.3350 ? ? 4.7717 ?-0.489 0.628274 ? ?

A ? ? ? ? ? ? 0.6113 ? ? 0.1476 ? 4.143 0.000271 ***

---

Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.831 on 29 degrees of freedom

Multiple R-squared: ?0.3718, Adjusted R-squared: ?0.3501?

F-statistic: 17.16 on 1 and 29 DF, ?p-value: 0.0002712

> #擬合直線

> plot(B~A,data=data)

> abline(pic)

> anova(pic)

Analysis of Variance Table

Response: B

? ? ? ? ? Df Sum Sq Mean Sq F value ? ?Pr(>F) ? ?

A ? ? ? ? ?1 251.85 251.846 ?17.162 0.0002712 ***

Residuals 29 425.57 ?14.675 ? ? ? ? ? ? ? ? ? ? ?

---

Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> ##擬合直線診斷-殘差分析:標(biāo)準(zhǔn)方法

> op<-par(mfrow=c(2,2))

> plot(pic)

> par(op)

> ###改進(jìn)的回歸擬合診斷方法：

> ##正態(tài)性檢驗(yàn)

> qqPlot(pic,id.method="identify",simlabels=row.names(data),ulate=TRUE,main="Normal Q-Q Plot")

[1] 14 19

> ###殘差獨(dú)立性檢驗(yàn) p>0.05,獨(dú)立

> library(car)

> durbinWatsonTest(pic)

?lag Autocorrelation D-W Statistic p-value

? ?1 ? ? ? 0.4027404 ? ? ?1.111413 ? 0.008

?Alternative hypothesis: rho != 0

> ##殘差線性?

> crPlots(pic)

> ##同方差性 >0.05,滿足方差不變

> ncvTest(pic)

Non-constant Variance Score Test?

Variance formula: ~ fitted.values?

Chisquare = 1.413038, Df = 1, p = 0.23455

> spreadLevelPlot(pic)

Suggested power transformation: ?3.43003?

> ##離群值

> outlierTest(pic)

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

? ?rstudent unadjusted p-value Bonferonni p

14 2.741071 ? ? ? ? ? 0.010545 ? ? ? 0.3269

> ##強(qiáng)影響點(diǎn)

> avPlots(pic)

> ####將離散點(diǎn)、杠桿值术裸、強(qiáng)影響點(diǎn)做到一個圖中 縱坐標(biāo)絕對值大于2為離群點(diǎn)倘是；水平大于0.2或者0.3的是高杠桿值；

> #圓圈大小與影響成正比

> influencePlot(pic,main="Influence Plot",sub="Circle size is proportional to Cook`s distance")

? ? ? ?StudRes ? ? ? ?Hat ? ? ? ?CookD

1 ?-1.69284861 0.10495836 1.578707e-01

7 ? 1.02479319 0.12721355 7.640411e-02

14 ?2.74107120 0.05599694 1.819728e-01

19 ?2.46795140 0.03819278 1.028721e-01

31 ?0.01449139 0.21178329 2.921952e-05

> ## 驗(yàn)證 可信度

> anova(pic)

Analysis of Variance Table

Response: B

? ? ? ? ? Df Sum Sq Mean Sq F value ? ?Pr(>F) ? ?

A ? ? ? ? ?1 251.85 251.846 ?17.162 0.0002712 ***

Residuals 29 425.57 ?14.675 ? ? ? ? ? ? ? ? ? ? ?

---

Signif. codes: ?0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#廣義線性模型-logistic 與泊松回歸 Y值類型為0/1或者一些計(jì)數(shù)型的

#logistics-AER包袭艺，婚外情數(shù)據(jù)affairs搀崭，各種因素是不是跟婚外情有關(guān)系

data(Affairs,package = "AER")

#變量x的統(tǒng)計(jì)分析

summary(Affairs)

#相應(yīng)變量y的統(tǒng)計(jì)數(shù)值

table(Affairs$affairs)

#這時(shí)候我們需要把婚外情變化成二值型數(shù)據(jù)0/1? 只關(guān)心有沒有? ynaffair可以任意取名字，新增的一列而已

Affairs$ynaffair[Affairs$affairs>0]<-1

Affairs$ynaffair[Affairs$affairs==0]<-0

Affairs$ynaffair<-factor(Affairs$ynaffair,levels = c(0,1),labels = c("No","Yes"))

table(Affairs$ynaffair)

##此時(shí)因?yàn)閅變?yōu)槎敌妥兞?/1匹表，可以用logistic回歸了

fit<-glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating,data=Affairs,family=binomial())

summary(fit)

#去除擬合不好的再擬合

fit_reduced<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())

summary(fit_reduced)

##因?yàn)閒it_reduced擬合是fit擬合的子集门坷，可以用卡方檢驗(yàn)的anova進(jìn)行比較? p>0.05，沒有差別袍镀，因此可以刪去無關(guān)的數(shù)據(jù)

anova(fit_reduced,fit,test = "Chisq")

##可以用coef單獨(dú)看看系數(shù)

coef(fit_reduced)

##方法一：用優(yōu)勢比解釋系數(shù)

##因?yàn)轫憫?yīng)變量是Y=1的log值默蚌，不好解釋，可以e-指數(shù)化后恢復(fù)原來的值苇羡；保持年齡什么的不變绸吸，婚齡增加一年，婚外情

? ? ##優(yōu)勢比乘于1.106.设江。锦茁。因此大于1，上升叉存，小于1下降優(yōu)勢比? 若分析幾年码俩，就乘于幾年就好了*n --n年

exp(coef(fit_reduced))

##系數(shù)置信區(qū)間

confint(fit_reduced)

##指數(shù)化的結(jié)果--一般用這個

exp(confint(fit_reduced))

##方法二：用概率解釋系數(shù)

##創(chuàng)建感興趣的數(shù)據(jù)集

testdata<-data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),

? ? ? ? ? ? ? ? ? ? religiousness=mean(Affairs$religiousness))

testdata

testdata$prob<-predict(fit_reduced,newdata = testdata,type = "response")

testdata

##過度離勢的判斷，避免不準(zhǔn)確的顯著性檢驗(yàn)

fit_reduced<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())

fit.od<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family = quasibinomial())

#檢驗(yàn) >0.34,沒有過度離勢歼捏，之前的檢驗(yàn)顯著可信稿存；

pchisq(summary(fit.od)$dispersion*fit_reduced$df.residual,fit_reduced$df.residual,lower=F)

##常見的回歸診斷可以用線性回歸診斷也可以用下面的

##一般方法：

op<-par(mfrow=c(2,2))

plot(fit_reduced)

par(op)

###本章的方法? 預(yù)測值與殘差值

plot(predict(fit_reduced,type = "response"),residuals(fit_reduced,type="deviance"))

###異常值--hat value? 學(xué)生化殘差，cook距離

plot(hatvalues(fit_reduced))

plot(rstudent(fit_reduced))

plot(cooks.distance(fit_reduced))

###泊松回歸分析---robust包

library(robust)

data(breslow.dat,package = "robust")

names(breslow.dat)

summary(breslow.dat[c(6,7,8,10)])

##基本圖形描述

op<-par(no.readonly = T)

par(mfrow=c(1,2))

attach(breslow.dat)

hist(sumY,breaks = 20,xlab = "Seizure Count",main = "Distribution of Seizures")

boxplot(sumY~Trt,xlab="Treatment",main="Group Comparisons")

par(op)

###回歸擬合：

fit_BS<-glm(sumY~Base+Age+Trt,family = poisson(),data=breslow.dat)

summary(fit_BS)

###系數(shù)分析

coef(fit_BS)

##指數(shù)化

exp(coef(fit_BS))

##過度離勢分析

##殘差偏差與殘差自由度比值? 遠(yuǎn)遠(yuǎn)大于1? 說明過度離勢

deviance(fit_BS)/df.residual(fit_BS)

#或用qcc包中進(jìn)行過度離勢分析? p<0.05 進(jìn)一步說明過度離勢

library(qcc)

qcc.overdispersion.test(breslow.dat$sumY,type="poisson")

##將poisson換成quasipoisson 重新擬合

fit<-glm(sumY~Base+Age+Trt,family = quasipoisson(),data=breslow.dat)

summary(fit)

##########################################################

##第五題? 回歸分析

data_5<-read.table("final5/diabetes.txt",header = T,sep = "\t")

data_5<-as.data.frame(data_5)

boxplot(data_5)

library(car)

scatterplotMatrix(data_5,spread=F,main="Scatter Plot Matrix")

##擬合

fit_5<-lm(Y~X1+X2+X3+X4,data = data_5)

summary(fit_5)

#回歸診斷--標(biāo)準(zhǔn)方法

op<-par(mfrow=c(2,2))

plot(fit_5)

par(op)

par(par(mfrow=c(2,2)))

##離群點(diǎn)

outlierTest(fit_5)

##逐步回歸分析 P大于0.05可以刪去x1 x2

fit_5_new<-lm(Y~X3+X4,data = data_5)

anova(fit_5,fit_5_new)

##交互作用

fit_5_ab<-lm(Y~X3+X4+X3:X4,data = data_5)

fit_5_ab

##交互結(jié)果展示

library(effects)

plot(effect("X3:X4",fit_5_ab),multiline=T)

##線性回歸包括幾個方面：數(shù)據(jù)觀察瞳秽，初步擬合瓣履，交互作用，多重共線性练俐，回歸診斷袖迎，擬合中出現(xiàn)的離群點(diǎn)，杠桿值腺晾，強(qiáng)影響點(diǎn)燕锥，刪除后，重新擬合悯蝉，模型比較--

##之后可能要預(yù)測归形，因此需要數(shù)據(jù)交叉驗(yàn)證分成2部分，一部分?jǐn)M合泉粉，一部分驗(yàn)證连霉。