数据集与问题描述
数据集测量了56条欧洲河鲈(Perca fluviatilis)的六项指标
- Weight Weight of the fish (in g)
- Length1 Length from the nose to the beginning of the tail (in cm)
- Length2 Length from the nose to the notch of the tail (in cm)
- Length3 Length from the nose to the end of the tail (in cm)
- Height Maximal height (in cm)
- Width Maximal width (in cm)
我们想知道鱼的重量和其他指标的关系。
数据集可以在 fish.data下载
相关库的加载和数据的读入
1
2
3
4
5
6
7
8
| requirements <- c('car', 'lars', 'tidyr', 'GGally')
for (package in requirements) {
if (!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}
fishData <- read.table('fish.data', header=T)
|
探索性数据分析
通常拿到一个问题首先要做的就是对数据的探索性分析。惨痛的教训告诉我们，这一步绝不应该是走流程。
1
2
3
4
5
6
7
8
9
10
11
12
13
| dim(fishData)
## [1] 56 6
head(fishData)
## Weight Length1 Length2 Length3 Height Width
## 1 5.9 7.5 8.4 8.8 2.1120 1.4080
## 2 32.0 12.5 13.7 14.7 3.5280 1.9992
## 3 40.0 13.8 15.0 16.0 3.8240 2.4320
## 4 51.5 15.0 16.2 17.2 4.5924 2.6316
## 5 70.0 15.7 17.4 18.5 4.5880 2.9415
## 6 100.0 16.2 18.0 19.2 5.2224 3.3216
ggpairs(fishData, pch=19, cex=.2) # see the picture below
n <- dim(fishData)[1] # 56
p <- dim(fishData)[2] # 6
|
从图中可以清晰的看到两点
- 五个自变量间高度相关，这和预期完全一致
- 数据出现了双峰的分布，仔细观察可以看到数据间存在一段gap(事实上，这是我们后期才发现的特性，因为当时只是使用了
pairs验证了自变量的高度相关性，漏掉了这一值得注意的信息)
尝试建模
一个显然不合理的模型
我们先尝试将所有变量拿来预测，直观上其结果必定会十分糟糕，我们只是验证一下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
| simpleModel <- lm(Weight~.,fishData)
summary(simpleModel)
## Call:
## lm(formula = Weight ~ ., data = fishData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -125.03 -52.44 -26.28 34.73 301.66
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -556.59 60.66 -9.175 2.68e-12 ***
## Length1 -3.13 57.88 -0.054 0.9571
## Length2 -38.50 88.55 -0.435 0.6656
## Length3 42.92 60.09 0.714 0.4784
## Height 65.66 30.00 2.189 0.0333 *
## Width 64.90 36.77 1.765 0.0836 .
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 87.11 on 50 degrees of freedom
## Multiple R-squared: 0.9429, Adjusted R-squared: 0.9372
## F-statistic: 165.2 on 5 and 50 DF, p-value: < 2.2e-16
vif(simpleModel)
## Length1 Length2 Length3 Height Width
## 1780.07795 4626.41486 2376.77263 54.03761 30.85782
|
可以从VIF看到自变量间高度相关（三个长度的VIF更是达到了令人发指的四位数），而五个系数的方差都很大，基本没有什么解释性，这就是之前谈到的多重共线性的危害。
使用Lasso选择变量
我们知道由于Lasso可以起到选择变量的作用，因为其切点更有可能出现在坐标轴上。我们使用lars实现这一功能
1
2
3
4
| plot(lars(fishData[,-1] %>% data.matrix(), fishData[,1]))
vif(lm(Weight~Height+Width, fishData))
## Height Width
## 29.56642 29.56642
|
从图中可以看到，随着惩罚的逐渐降低，变量4和变量5依次加入到模型中，我们检查这样的模型的VIF，发现仍然高达30。因此即使是双变量都存在着严重的共线性问题，这从根本上断绝了多元回归的路子。
PCA降维
既然现在变量间高度相关，使用PCA提取出一个“综合体型”这样的想法是十分自然的。我们先使用PCA看下效果
1
2
| pca <- prcomp(fishData[,-1], center=F, scale=F)
screeplot(pca, type='lines')
|
毫不意外的，第一主成分占走了绝大多数方差，因此我们只能选用第一主成分。然后我们通过散点图看下主成分和因变量之间的关系。
1
2
3
| y <- fishData[,1]
x <- -pca$x[,1]
plot(x, y, pch=19, cex=.2)
|
显然线性模型是不合适的，从直观上体重也不会随身高线性增长，相反他们应该有一个幂的关系。因此我们考虑对双放取对数后转化为线性模型来求解。
但是，现实中使用PCA是一个很麻烦的问题，你还得手动测量五个长度，然后根据一堆奇奇怪怪的系数（这个系数也不见得稳定，因为样本bias）换算得出。注意到其和Length是高度相关的，因此不如直接用某个长度更为方便直接。三个Length用哪个呢？个人更倾向于使用Length2，因为其测量相对简单：量头部到尾巴缺口的距离即可，不像尾巴的头部或尾部难以定位。因此后面我们都只使用Length2代替x。
对数化模型
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
| plot(y=log(y), x=log(x), pch=19, cex=.2)
log.fit <- lm(log(y)~log(x))
summary(log.fit)
## Call:
## lm(formula = log(y) ~ log(x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21311 -0.08150 -0.00789 0.04987 0.38451
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.86097 0.15158 -32.07 <2e-16 ***
## log(x) 3.15296 0.04606 68.46 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 0.1177 on 54 degrees of freedom
## Multiple R-squared: 0.9886, Adjusted R-squared: 0.9884
## F-statistic: 4686 on 1 and 54 DF, p-value: < 2.2e-16
|
取对数后可以看到一个相当漂亮的线性关系，其拟合效果也相当不错，\(R^2\)达到了0.9886。事实上这个模型基本上可以收工了。但当观察到\(\log(x)\)的系数为\(3.15\)时，我们不禁有一个大胆的想法——如果用简单的三次方模型是否可行呢？
当然，冷静看一下标准差，不难看出\(3\)不太可能落在\(95\%\)的置信区间里。让我们用代码尝试一下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| summary(lm(log(y)-3*log(x)~log(x)))
## Call:
## lm(formula = log(y) - 3 * log(x) ~ log(x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21311 -0.08150 -0.00789 0.04987 0.38451
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.86097 0.15158 -32.068 < 2e-16 ***
## log(x) 0.15296 0.04606 3.321 0.00161 **
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 0.1177 on 54 degrees of freedom
## Multiple R-squared: 0.1696, Adjusted R-squared: 0.1542
## F-statistic: 11.03 on 1 and 54 DF, p-value: 0.001613
|
\(\beta=3\)的p-value为0.00161，算是不太可能了。当然我们也可以argue说因为数据量太小可能有bias，因此我们就使用一个三次方模型使得模型更好解释，这就见仁见智了。
尝试更高阶的模型
之前我们得到了一个\(3.12\)阶的模型，然后发现降到\(3\)阶的想法不太被支持，但我们确实更喜欢整数次方的拟合结果，那么\(y=a_0+a_1x+a_2x^2+a_3x^3\)怎么样呢？
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
| x2 <- x^2
x3 <- x^3
summary(lm(y~x3+x2))
## Call:
## lm(formula = y ~ x3 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -154.792 -28.543 -3.992 18.650 249.137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -60.410113 38.040283 -1.588 0.11822
## x3 0.009296 0.002942 3.160 0.00261 **
## x2 0.206276 0.139342 1.480 0.14470
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 60.65 on 53 degrees of freedom
## Multiple R-squared: 0.9707, Adjusted R-squared: 0.9696
## F-statistic: 876.9 on 2 and 53 DF, p-value: < 2.2e-16
summary(lm(y~x3+x))
## Call:
## lm(formula = y ~ x3 + x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.74 -24.75 -4.65 18.25 254.74
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -76.144274 65.819272 -1.157 0.253
## x3 0.012291 0.001292 9.514 4.66e-13 ***
## x 3.840499 3.597097 1.068 0.291
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 61.24 on 53 degrees of freedom
## Multiple R-squared: 0.9701, Adjusted R-squared: 0.969
## F-statistic: 859.7 on 2 and 53 DF, p-value: < 2.2e-16
|
可以看到当\(x^3\)引入后，很难继续引入\(x\)与\(x^2\)。当然，其实这么直接粗暴的引入\(x\)和\(x^2\)是有问题的——会带来共线性问题。但这里我们只是做一些探索性的尝试，就不去标准化了。
数据标准化
当然，尝试下对数据标准化也不是不行，只是在这个数据上结果并没有显著改善，而且解释性也会变得很差（试想你要如何解释标准化的系数？只能说我们有一个“标准”的鱼，然后根据体型上对标准鱼的偏差计算出重量上对标准鱼的偏差，再换算会鱼的重量。因此这里只放相关的代码，有兴趣的读者自行尝试。
1
2
3
4
5
6
7
8
| nx1 <- scale(x)
nx2 <- nx1^2
nx3 <- nx1^3
ny <- scale(y)
plot(nx1, ny, pch=19, cex=.2)
nfit <- lm(ny~nx1)
step(nfit, scope=list(lower=ny~1, upper=ny~nx1+nx2+nx3))
summary(lm(ny~nx1+nx2))
|
阶段性成果
经过思考，我们最后还是无法抵挡整数次方带来的诱惑，因此采取\(y=a+bx^3\)这样的模型。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
| summary(lm(y~x3))
## Call:
## lm(formula = y ~ x3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -153.502 -17.969 -5.177 9.488 262.964
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.1415073 12.4741568 -0.573 0.569
## x3 0.0136251 0.0003291 41.398 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 61.32 on 54 degrees of freedom
## Multiple R-squared: 0.9695, Adjusted R-squared: 0.9689
## F-statistic: 1714 on 1 and 54 DF, p-value: < 2.2e-16
|
注意到，截距项的p-value达到了0.693，这意味着我们可以抛弃掉这一参数而得到更精炼的模型
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
| newModel <- lm(y~x3+0)
summary(newModel)
## Call:
## lm(formula = y ~ x3 + 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -149.801 -23.571 -11.324 9.543 261.914
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x3 0.0134831 0.0002149 62.75 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 60.94 on 55 degrees of freedom
## Multiple R-squared: 0.9862, Adjusted R-squared: 0.986
## F-statistic: 3938 on 1 and 55 DF, p-value: < 2.2e-16
plot(x3, y, xlab="x^3", pch=19, cex=.2)
lines(abline(newModel))
|
值得注意的时，使用\(\log(y)=3*\log(x)+b\)也能得到一个模型，这模型和上面直接三次方你和出来的模型是不一样的。选用哪种模型也是一个见人见智的事情。直接\(y=ax^3\)的\(\mathrm{SSE}\)肯定更小，但这样拟合的模型自变量相对不太好——如果\(x\)相对还比较均匀，那么\(x^3\)到后期可能就相当稀疏了，其误差带来的影响也会被放大，因此此时我们可能更喜欢对数化后的线性模型。
1
2
3
4
| mean(log(y)-3*log(x))
# [1] -4.360281
exp(mean(log(y)-3*log(x)))
# [1] 0.0127748
|
得到的模型为 \(y=0.01277 x^3\)。
在下面的讨论中，我们都使用对数化的模型做进一步的改进。当然读者也可以尝试用\(y=ax^3\)直接拟合然后尝试下面的改进。
还可以更好
上述这个模型基本上没有问题了，简洁明了有针对性，物理解释也足够好，过零点简直是一个意外的bonus。但当做模型诊断的时候，我们发现了这样的问题。
1
2
3
4
5
6
| par(mfrow=(c(2,2)))
plot(x3, newModel$residuals, xlab="x^3", ylab="e", pch=19, cex=.2)
plot(predict(newModel), newModel$residuals, xlab="yhat", ylab="e", pch=19, cex=.2)
plot(newModel$residuals[-n], newModel$residuals[-1], xlab="e_i", ylab="e_i+1", pch=19, cex=.2)
qqnorm(newModel$residuals, pch=19, cex=.2)
sum(newModel$residuals^2)/(n-1)
|
可以看到一个是误差的正态性不太好，有重尾的趋势。但事实上正态性假设在现实中许多时候是可有可无的。
另一个问题就是随着\(X\)或者\(\widehat{Y}\)的增大，误差也会随之增大，但这事实上是因为\(x^3\)放大了这一效应，在对数视角上看，其误差并没有明显的增大，还是可以接受的。
但我们能不能做一些改进呢？还记得最开始做探索性分析的时候我们发现的现象么？我们可以做这样的推想：这里面的鱼有幼年期和成熟期之分（或者是两个亚种？），因此体型有较大的差异。这可能导致他们三次方的系数不一样。因此，我们可以使用一个因子isBig标志其是否为大鱼，然后做一个协方差分析。
经过观察，数据的gap在Height比较大，因此我们选择在此维度作为划分，且取了gap中一个中庸的值9.0作为划分。此时，我们设想的模型是
\[Y=(a+b\cdot\mathrm{isBig})x^3\]
你可以使用两种方法得到相关的参数，一个简单的方法自然是分别回归两个模型，然后拼凑出对应的系数，即
1
2
3
4
5
6
7
| isBig <- fishData$Height > 9
small.coef <- mean(log(y[!isBig])-3*log(x[!isBig]))
big.coef <- mean(log(y[isBig])-3*log(x[isBig]))
exp(small.coef)
# [1] 0.01225077
exp(big.coef)
# [1] 0.01386048
|
当然另一个做法是直接使用回归的方式，这样做的好处是可以让R帮我们算一下\(t\)检验等。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| final.model <- lm(log(y)-3*log(x)~isBig)
summary(final.model)
## Call:
## lm(formula = log(y) - 3 * log(x) ~ isBig)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.23500 -0.06384 -0.00319 0.04192 0.34607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.40217 0.01884 -233.716 < 2e-16 ***
## isBigTRUE 0.12345 0.03234 3.818 0.000349 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 0.1146 on 54 degrees of freedom
## Multiple R-squared: 0.2125, Adjusted R-squared: 0.198
## F-statistic: 14.58 on 1 and 54 DF, p-value: 0.0003491
|
可以看到这个因子的引入还是有用的。
至此，我们得到了我们使用的最终模型
\[\log(\mathrm{Weight})=(-4.40217+0.12345\mathrm{isBig}) + 3\log(\mathrm{Length2})\] \[\mathrm{Weight}=(0.0125+0.00161\mathrm{isBig}) (\mathrm{Length2})W^3\]
模型诊断
由于我们已经写死了三次幂的关系，因此一般的像检测\(X,\,Y\)的离群值意义相对不大，因此只做一些基础性的诊断分析。
1
| qqnorm(final.model$resuiduals)
|
至此，一个线性模型的分析任务就完成了，谢谢各位读者。