第九章 雙變量回歸與相關(guān)
x2yline
r Sys.Date()
知識(shí)清單
- 直線(xiàn)回歸
- 相關(guān)概念
- 求法
- 統(tǒng)計(jì)推斷
- 區(qū)間估計(jì)
- 直線(xiàn)相關(guān)
- 相關(guān)概念
- 求法
- 統(tǒng)計(jì)推斷
- 決定系數(shù)
- 直線(xiàn)回歸應(yīng)用注意事項(xiàng)
- 秩相關(guān)
- 兩回歸直線(xiàn)的比較
- 曲線(xiàn)擬合
1. 直線(xiàn)回歸
1.1 基本概念
- 用途
- 不是前幾章的單變量堵腹，而是兩個(gè)變量的關(guān)系
- 數(shù)值變量或有序分類(lèi)變量（秩相關(guān)）
- 直線(xiàn)方程及各值含義
- 方程 1
[
\hat{y}=a=bX
] - 方程 2
[
\mu_{Y|X}=\alpha+\beta X
]- $\hat{y}$表示X對(duì)應(yīng)Y的總體均數(shù)$\mu_{Y|X}$的一個(gè)樣本估計(jì)，稱(chēng)預(yù)測(cè)值
- a泉孩、b分別為$\alpha$和$\beta$的樣本估計(jì)留储，a稱(chēng)為常數(shù)項(xiàng)夹囚，b稱(chēng)回歸系數(shù)
- b的統(tǒng)計(jì)意義是當(dāng)X變化一個(gè)單位時(shí)Y的平均改變的估計(jì)值
- 方程 1
1.2 計(jì)算方法
1.2.1 理論計(jì)算
根據(jù)實(shí)測(cè)值Y與假定回歸直線(xiàn)上的估計(jì)值$\hat{y}$的縱向距離$Y-\hat{y}$（殘差）的平方和最小即最小二乘法可以推出
[
b=\frac{l_{XY}}{l_{XX}}=\frac{\sum{}{}(X-\bar{X})(Y-\bar{Y})}{\sum{}{}(X-\bar{X})^{2})}
]
[
a=\bar{Y}-b\bar{X}
]
式中$l_{XY}$為X與Y的離均差交叉乘積和，簡(jiǎn)稱(chēng)離均差積和
1.2.2 R語(yǔ)言計(jì)算
data9_1 <- haven::read_sav(
file="E:\\醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）\\各章例題SPSS數(shù)據(jù)文件\\例09-01.sav")
colnames(data9_1) <- c("age", "conc")
head(data9_1)
## # A tibble: 6 x 2
## age conc
## <dbl> <dbl>
## 1 13 3.54
## 2 11 3.01
## 3 9 3.09
## 4 6 2.48
## 5 8 2.56
## 6 10 3.36
line.model <- lm(conc~age, data=data9_1)
print(line.model)
##
## Call:
## lm(formula = conc ~ age, data = data9_1)
##
## Coefficients:
## (Intercept) age
## 1.6617 0.1392
line.model_summary <- summary(line.model)
1.3 統(tǒng)計(jì)推斷
H0：$\beta=0$，即兩變量直接無(wú)直線(xiàn)關(guān)系
H1：$\beta \neq 0$争群，即兩變量之間有線(xiàn)性關(guān)系
1.3.1 方差分析
公式
把Y的離均差平方和進(jìn)行分解峰弹，分為回歸平方和與殘差平方和店量，其自由度分別為1，n-2鞠呈。
[
Y-\bar{Y}=Y-\hat{Y}+\hat{Y}-\bar{Y}
]
[
\sum{}{}(Y-\bar{Y}){2}=\sum{}{}(\hat{y}-\bar{Y}){2} + \sum{}{}(Y-\hat{y})^{2}
]
[
SS_{總}(cāng)=SS_{回}+SS_{殘}
]
F值計(jì)算公式為（單側(cè)F檢驗(yàn)）：
[
F=\frac{SS_{回}/\nu_{回}}{SS_{殘}/\nu_{殘}},\ \nu_{回}=1,\ \nu_{殘}=n-2
]
$SS_{回}$計(jì)算公式可以化簡(jiǎn)為：
[
SS_{回}=bl_{XY}=\frac{l_{XY}{2}}{l_{XX}}=b{2}l_{XX}
]
R語(yǔ)言實(shí)現(xiàn)
- 手動(dòng)計(jì)算
# 殘差平方和
ss2 <- sum(line.model$residuals^2)
# 離均差平方和
ss0 <- var(data9_1$conc)*(nrow(data9_1)-1)
# 回歸平方和
ss1 <- ss0-ss2
# F統(tǒng)計(jì)量
f.statistic <- (ss1/1)/(ss2/(nrow(data9_1)-2))
# p值
p <- pf(f.statistic, lower.tail=FALSE, df1=1, df2=nrow(data9_1)-2)
cat("F statistic is ", f.statistic, "\np value is ", p, sep="")
## F statistic is 20.96842
## p value is 0.003773985
- 直接調(diào)用查看summary.lm對(duì)象里的f值并轉(zhuǎn)為p值
f_df1_df2 <- summary(line.model)$fstatistic
p_value <- pf(f_df1_df2[1], df1=f_df1_df2[2], df2=f_df1_df2[3], lower.tail=F)
cat(cat("F statistic is ", f_df1_df2[1], "\np value is ", p_value, sep=""))
## F statistic is 20.96842
## p value is 0.003773985
- 直接打印sammary.lm對(duì)象融师，最后一行信息即為其F值和對(duì)應(yīng)p值
summary(lm(conc~age, data=data9_1))
##
## Call:
## lm(formula = conc ~ age, data = data9_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21500 -0.15937 -0.00125 0.09583 0.30667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.66167 0.29700 5.595 0.00139 **
## age 0.13917 0.03039 4.579 0.00377 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.197 on 6 degrees of freedom
## Multiple R-squared: 0.7775, Adjusted R-squared: 0.7404
## F-statistic: 20.97 on 1 and 6 DF, p-value: 0.003774
1.3.2 t檢驗(yàn)
t值計(jì)算公式：
[
t=\frac{b-0}{S_},\ \nu=n-2
]
[
S_蚁吝=\frac{S_{Y\cdot X}}{\sqrt{l_{XX}}}
]
[
S_{Y\cdot X}=\sqrt{\frac{SS_{殘}}{n-2}}
]
$S_{Y\cdot X}$為回歸的剩余標(biāo)準(zhǔn)差旱爆，化簡(jiǎn)后有：
[
\sqrt{F}=t
]
R語(yǔ)言實(shí)現(xiàn)（雙側(cè)t檢驗(yàn)）：
lxx <- var(data9_1$age)*(nrow(data9_1)-1)
Sb <- sqrt(sum(line.model$residuals^2)/(nrow(data9_1)-2))/(sqrt(lxx))
t_statistic <- (line.model$coefficients[2]-0)/Sb
cat("t statistic is ", t_statistic, "\np value is ",
2*pt(t_statistic, df=nrow(data9_1)-2, lower.tail=FALSE),
sep="")
## t statistic is 4.579129
## p value is 0.003773985
1.3.3 區(qū)間估計(jì)
- 總體回歸系數(shù)$\beta$的可信區(qū)間
表示Y0的均數(shù)95%的置信區(qū)間
結(jié)合上述t檢驗(yàn)的公式，$\beta$的$1-\alpha$可信區(qū)間為：
[
b\pm t_{\alpha/2,\ \nu}\cdot S_窘茁
]
- 總體均數(shù)$\mu$的可信區(qū)間
$\hat{Y0}$會(huì)因樣本（擬合的曲線(xiàn)）而異怀伦，其抽樣誤差大小的標(biāo)準(zhǔn)誤：
[
S_{\hat{Y_{0}}}=S_{Y\cdot X}\sqrt{\frac{1}{n}+\frac{(X_{0}-\bar{X}){2}}{\sum{}{}(X-\bar{X}){2}}}
]
$\mu_{Y|X0}$的置信區(qū)間為：
[
\hat{Y_{0}} \pm t_{\alpha/2,\ \nu}\cdot S_{\hat{Y_{0}}}
]
R語(yǔ)言實(shí)現(xiàn)：
lxx <- var(data9_1$age)*(nrow(data9_1)-1)
# 剩余標(biāo)準(zhǔn)差
syx <- sqrt(sum(line.model$residuals^2)/(nrow(data9_1)-2))
Sb <- syx/(sqrt(lxx))
t_statistic <- (line.model$coefficients[2]-0)/Sb
Sy01 <- syx*sqrt(1/nrow(data9_1)+
(data9_1$age-mean(data9_1$age))^2/
(var(data9_1$age)*nrow(data9_1)-1))
# geom_smooth自動(dòng)加上了標(biāo)準(zhǔn)偏差即se=TRUE
library(ggplot2)
p <- ggplot(data9_1, aes(age, conc))+
xlab("age")+ ylab("concentration")+
geom_point(size=1.5)+
geom_smooth(method="lm", se=TRUE)+
theme_classic()
print(p)
p <- p + geom_smooth(aes(x=age, y=Sy01+
age*line.model$coefficients[2]+
line.model$coefficients[1]), color="red",
se=FALSE, method="loess")+
geom_smooth(aes(x=age, y=-Sy01+
age*line.model$coefficients[2]+
line.model$coefficients[1]), color="red",
se=FALSE, method="loess")
print(p)
- 個(gè)體Y值的預(yù)測(cè)區(qū)間
表示Y0值的95%置信區(qū)間范圍
[
S_{Y_{0}}=S_{Y\cdot X}\sqrt{1+\frac{1}{n}+\frac{(X_{0}-\bar{X}){2}}{\sum{}{}(X-\bar{X}){2}}}
]
R語(yǔ)言實(shí)現(xiàn)：
lxx <- var(data9_1$age)*(nrow(data9_1)-1)
# 剩余標(biāo)準(zhǔn)差
syx <- sqrt(sum(line.model$residuals^2)/(nrow(data9_1)-2))
Sb <- syx/(sqrt(lxx))
t_statistic <- (line.model$coefficients[2]-0)/Sb
Sy02 <- syx*sqrt(1+1/nrow(data9_1)+
(data9_1$age-mean(data9_1$age))^2/
(var(data9_1$age)*nrow(data9_1)-1))
p + geom_smooth(aes(x=age, y=Sy02+
age*line.model$coefficients[2]+
line.model$coefficients[1]), color="green",
se=FALSE, method="loess")+
geom_smooth(aes(x=age, y=-Sy02+
age*line.model$coefficients[2]+
line.model$coefficients[1]), color="green",
se=FALSE, method="loess")+
geom_vline(xintercept=mean(data9_1$age), lwd=1, color="yellow", linetype=2)+
annotate("text", x=mean(data9_1$age), y=3.5, label="X_bar")+
guides(color = guide_legend(title = "LEFT", title.position = "left"))
2. 直線(xiàn)相關(guān)
2.1 相關(guān)概念
- 又稱(chēng)簡(jiǎn)單相關(guān)，用于雙變量正態(tài)分布
- 相關(guān)系數(shù)（coefficient of correlation）又稱(chēng)Pearson積差相關(guān)系數(shù)（coefficient of product-moment correlation）山林，符號(hào)$r$代表樣本相關(guān)系數(shù)房待，符號(hào)$p$代表總體相關(guān)系數(shù)。
2.2 計(jì)算公式
以符號(hào)$r$表示樣本相關(guān)系數(shù)驼抹，符號(hào)$\rho$表示總體相關(guān)系數(shù)桑孩，$r$是$rho$的估計(jì)，與b不同框冀，它沒(méi)有單位
[
r=\frac{\sum{}{}(X-\bar{X})(Y-\bar{Y})}{\sqrt{\sum{}{}(X-\bar(X)){2}}\sqrt{\sum{}{}(Y-\bar{Y}){2}}}=\frac{l_{XY}}{\sqrt{l_{XY}l_{YY}}}
]
R語(yǔ)言實(shí)現(xiàn)：
data9_5 <- haven::read_sav(
file="E:\\醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）\\各章例題SPSS數(shù)據(jù)文件\\例09-05.sav")
colnames(data9_5) <- c("number", "age", "v")
# 公式法
r <- sum((data9_5$age-mean(data9_5$age))*(data9_5$v-mean(data9_5$v)))/
sqrt(var(data9_5$age)*(nrow(data9_5)-1)*var(data9_5$v)*(nrow(data9_5)-1))
print(r)
## [1] 0.8754315
# 包法
cor(data9_5$age, data9_5$v)
## [1] 0.8754315
2.3 統(tǒng)計(jì)推斷
2.3.1 t檢驗(yàn)
H0：$\rho=0$
[
S_{r}=\sqrt{\frac{1-r^{2}}{n-2}},\ \nu=n-2
]
[
t=\frac{r-0}{S_{r}}
]
R語(yǔ)言實(shí)現(xiàn)：
# 方法1：
Sr <- sqrt((1-r^2)/(nrow(data9_5)-2))
t_statistic <- (r-0)/Sr
pt(t_statistic, lower.tail=FALSE, df=nrow(data9_5)-2)*2
## [1] 1.910939e-05
# 方法2：
cor.test(data9_5$v, data9_5$age)
##
## Pearson's product-moment correlation
##
## data: data9_5$v and data9_5$age
## t = 6.5304, df = 13, p-value = 1.911e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6584522 0.9580540
## sample estimates:
## cor
## 0.8754315
2.3.2 可信區(qū)間
由于相關(guān)系數(shù)的抽樣分布在$\rho$不等于0的情況下呈偏態(tài)分布流椒，所以不能用t分布直接計(jì)算，向進(jìn)行變量變換明也，使其服從正態(tài)分布在計(jì)算可信區(qū)間
對(duì)r作z反雙曲正切函數(shù)變換：
[
z=tanh^{-1}r或z=\frac{1}{2}ln\frac{1+r}{1-r}
]近似計(jì)算z的可信區(qū)間
[(z-u_{\alpha/2/\sqrt{n-3}},\ z+u_{\alpha/2/\sqrt{n-3}})
]變換z為r
[
r=tanhz或r=\frac{e{2z}-1}{e{2z}+1}
]
R語(yǔ)言實(shí)現(xiàn)
# 方法1：
z <- atanh(r)
u_0.05_2 <- qnorm(0.975, mean=0, sd=1)
r1 <- tanh(z-u_0.05_2/sqrt(nrow(data9_5)-3))
r2 <- tanh(z+u_0.05_2/sqrt(nrow(data9_5)-3))
cat("CL is ", r1, "~", r2, sep="")
## CL is 0.6584522~0.958054
# 方法2：
cor.test(data9_5$age, data9_5$v, conf.level=0.95, alternative="two.sided")
##
## Pearson's product-moment correlation
##
## data: data9_5$age and data9_5$v
## t = 6.5304, df = 13, p-value = 1.911e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6584522 0.9580540
## sample estimates:
## cor
## 0.8754315
2.4 決定系數(shù)
定義：回歸平方和與總平方和之比宣虾，計(jì)算公式為
[
R{2}=\frac{SS_{回}}{SS_{總}(cāng)}=\frac{l_{XY}{2}/l_{XX}}{l_{YY}}
]
對(duì)于雙變量回歸分析，$R{2}$即$r{2}$温数，處可以概括擬合效果外绣硝，還可以作假設(shè)檢驗(yàn)
[
F=\frac{R{2}}{(1-R{2})(n-2)}=\frac{SS_{回}}{SS_{殘}/\nu_{殘}}=\frac{MS_{回}}{MS_{殘}}
]
R語(yǔ)言實(shí)現(xiàn)：
data9_5 <- haven::read_sav(
file="E:\\醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）\\各章例題SPSS數(shù)據(jù)文件\\例09-05.sav")
colnames(data9_5) <- c("number", "age", "v")
line.model <- lm(v~age, data=data9_5)
r.squared <- summary(line.model)$r.squared
print(r.squared)
## [1] 0.7663803
adj.r.squared <- summary(line.model)$adj.r.squared
print(adj.r.squared)
## [1] 0.7484095
r.squared <- 1-var(line.model$residuals)*(nrow(data9_5)-1)/
(var(data9_5$v)*(nrow(data9_5)-1))
print(r.squared)
## [1] 0.7663803
f.statistic <- (r.squared)/((1-r.squared)/(nrow(data9_5)-2))
print(f.statistic)
## [1] 42.64598
summary(line.model)$fstatistic
## value numdf dendf
## 42.64598 1.00000 13.00000
3. 直線(xiàn)回歸應(yīng)用注意事項(xiàng)
根據(jù)目的選擇變量及統(tǒng)計(jì)方法（自變量和因變量，如重測(cè)信度評(píng)價(jià)的相關(guān)系數(shù)帆吻，r應(yīng)達(dá)到0.40以上
進(jìn)行相關(guān)域那，回歸分析前應(yīng)繪制散點(diǎn)圖，離群值
-
結(jié)果解釋
- 相關(guān)系數(shù)或回歸系數(shù)的絕對(duì)值反映密切程度
- p值越小越有理由認(rèn)為變量間的直線(xiàn)關(guān)系存在
-
殘差圖觀(guān)察是否符合模型假設(shè)的條件：自變量與因變量關(guān)系為線(xiàn)性，誤差服從均數(shù)為0的正態(tài)分布次员，且方差相等败许，各觀(guān)測(cè)獨(dú)立(殘差圖橫坐標(biāo)為$\hat{Y}$或者X，縱坐標(biāo)為殘差)
- 正常為在y=0處對(duì)稱(chēng)分布淑蔚，并且左右對(duì)稱(chēng)
- 離群值會(huì)與群體遠(yuǎn)離
- 喇叭狀（左右不對(duì)稱(chēng)）說(shuō)明方差不齊（須穩(wěn)定化處理）
- 呈曲線(xiàn)則可能是符合曲線(xiàn)模型
- 呈直線(xiàn)說(shuō)明殘差與時(shí)間存在相關(guān)
R語(yǔ)言實(shí)現(xiàn)殘差圖
data9_5 <- haven::read_sav(
file="E:\\醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）\\各章例題SPSS數(shù)據(jù)文件\\例09-05.sav")
colnames(data9_5) <- c("number", "age", "v")
line.model <- lm(v~age, data=data9_5)
# +geom_point(aes(age, v), size=1.5)+
# geom_smooth(aes(age, v), method="lm")+
library(ggplot2)
ggplot(data9_5)+
geom_point(aes(x=predict(line.model), y=residuals(line.model)), color="red", size=2)+
ylab("residuals")+xlab(latex2exp::TeX("$\\hat{Y}$"))
殘差圖沒(méi)有明顯的偏倚趨勢(shì)（各區(qū)域殘差的變異程度大致相同）市殷，說(shuō)明殘差至少在一定的范圍內(nèi)是恒定的，該線(xiàn)性模型的效果基本還行
4. 秩相關(guān)
4.1 適用條件
不服從雙變量正態(tài)分布刹衫，而不宜做積差相關(guān)分析（散點(diǎn)圖或統(tǒng)計(jì)表看出）
總體分布類(lèi)型未知
原始數(shù)據(jù)用等級(jí)表示
4.2 Spearman秩相關(guān)
等級(jí)相關(guān)系數(shù)公式
$d$是指兩個(gè)變量的秩差醋寝，完全正相關(guān)則$\sum{}{}d_{i}{2}$有最小值為0，完全負(fù)相關(guān)則$\sum{}{}d_{i}{2}$有最大值為$\frac{n(n{2}-1}{3}$带迟，0相關(guān)則則$\sum{}{}d_{i}{2}=\frac{0+\frac{n(n{2}-1}{3}}{2}=\frac{n(n{2}-1}{6}$音羞，有以下公式
[
r_{s}=1-\frac{6\sum{}{}d{2}}{n(n{2}-1)}
]
$r_{s}$介于-1與1之間，負(fù)數(shù)則為負(fù)相關(guān)仓犬，正數(shù)則為正相關(guān)嗅绰，0則0相關(guān)
統(tǒng)計(jì)推斷
樣本等級(jí)相關(guān)系數(shù)$r_{s}$是總體相關(guān)系數(shù)$\rho_{s}$的估計(jì)值，檢驗(yàn)$\rho_{s}$是否不為0可以用查表法搀继，如果n大于50窘面，可以用u檢驗(yàn)，其中$u=r_{s}\sqrt{n-1}$叽躯，查u界值確定p
相同秩較多的情況下财边，需要校正，也可以不校正直接進(jìn)行秩的pearson相關(guān)系數(shù)的計(jì)算
R語(yǔ)言實(shí)現(xiàn)
data9_8 <- haven::read_sav(
file="E:\\醫(yī)學(xué)統(tǒng)計(jì)學(xué)（第4版）\\各章例題SPSS數(shù)據(jù)文件\\例09-08.sav")
colnames(data9_8) <- c("number", "X", "Y")
# 默認(rèn)method為pearson点骑，如果有相同秩酣难，會(huì)自動(dòng)校正
1-sum((rank(data9_8$X)-rank(data9_8$Y))^2)*6/(nrow(data9_8)^3-nrow(data9_8))
## [1] 0.9050568
cor(data9_8$X, data9_8$Y, method="spearman")
## [1] 0.9050568
cor.test(data9_8$X, data9_8$Y, method="spearman")
##
## Spearman's rank correlation rho
##
## data: data9_8$X and data9_8$Y
## S = 92, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.9050568
# 模擬有相同秩的情況
a <- c(1, 2, 3.5, 3.5, 5, 6)
b <- c(2, 1, 3, 4, 5.5, 5.5)
1-sum((rank(a)-rank(b))^2)*6/(length(a)^3-length(a))
## [1] 0.9142857
cor(a, b, method="spearman")
## [1] 0.9117647
cor函數(shù)還有一個(gè)use參數(shù)來(lái)處理數(shù)據(jù)缺失的情況（NA)，默認(rèn)為use="all.obs"黑滴，該情況下鲸鹦，如果數(shù)據(jù)有缺失，則報(bào)錯(cuò)跷跪，可以修改為use="complete.obs"，把含缺失數(shù)據(jù)的那一列刪除后再運(yùn)行齐板，或者修改為use="pairwise.complete.obs"吵瞻，把含缺失數(shù)據(jù)的那一行刪除后再運(yùn)行
