我有兴趣在 GAMM 中单独测试线性趋势和非线性分量。
我遵循了 mgcv 文档中的示例以及 stackexchange 上相关问题的答案,以使用
m=c(2,0)
和默认的 tp
平滑。该示例符合我的预期:(1) x 的参数系数估计的斜率非常接近简单线性回归的估计; (2) 平滑看起来去趋势且相关(非常小)大小。
但是,当我尝试用我的数据执行此操作时,结果对我来说毫无意义。看起来线性分量是任意的(有时甚至没有正确的符号,如下例所示)并且平滑没有去趋势(它显然包含强线性分量)。
我想知道我是否误解了这应该做什么,或者我是否以某种方式错误地编码了。由于文档示例工作正常,我不认为这是 mgcv 中的错误。
我还能如何在 GAMM 中独立于非线性分量来测试线性趋势?
这是 mgcv 文档中的相关示例:
require(mgcv)
n <- 100
set.seed(2)
x <- runif(n)
y <- x + x^2*.2 + rnorm(n) *.1
ed <- data.frame(y,x)
mg <- gam(y~s(x,m=c(2,0))+x,data=ed,method="REML")
据我了解,这从平滑中删除了线性分量,以便可以单独估计。
这是 gam 和 lm 的斜率估计值,非常接近:
> summary(mg)
Family: gaussian
Link function: identity
Formula:
y ~ s(x, m = c(2, 0)) + x
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02133 0.05315 -0.401 0.689
x 1.18249 0.10564 11.193 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(x) 0.9334 8 0.304 0.076 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.91 Deviance explained = 91.1%
-REML = -70.567 Scale est. = 0.012767 n = 100
> summary(lm(y~x,ed))
Call:
lm(formula = y ~ x, data = ed)
Residuals:
Min 1Q Median 3Q Max
-0.24082 -0.07942 -0.01136 0.09430 0.22836
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03019 0.02212 -1.365 0.175
x 1.20052 0.03851 31.175 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1144 on 98 degrees of freedom
Multiple R-squared: 0.9084, Adjusted R-squared: 0.9075
F-statistic: 971.9 on 1 and 98 DF, p-value: < 2.2e-16
这是带有 lm(红色)和 gam(蓝色)的叠加预测的数据图:
plot(y~x,ed)
abline(lm(y~x,ed),col="red")
points(ed$x,predict(mg,ed),col="blue")
这是平滑的,它看起来完全去趋势化并且非常小,与图表一致:
plot.gam(mg)
这是使用我的数据进行的相同分析。
read.table("ev.txt") -> ev
me <- gam(logart~s(item,m=c(2,0))+item,data=ev,method="REML")
由 gam 估计的项目的线性斜率是负值 (-0.11),尽管实际斜率是正值 (+0.32):
> summary(me)
Family: gaussian
Link function: identity
Formula:
logart ~ s(item, m = c(2, 0)) + item
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.305740 0.003587 1758.100 <2e-16 ***
item -0.112921 0.049174 -2.296 0.0217 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(item) 6.125 8 3.583 3.15e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.0298 Deviance explained = 3.14%
-REML = -842.54 Scale est. = 0.039768 n = 4453
> summary(lm(logart~item,ev))
Call:
lm(formula = logart ~ item, data = ev)
Residuals:
Min 1Q Median 3Q Max
-0.65009 -0.12830 0.00627 0.13301 0.61885
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.311575 0.003001 2103.47 <2e-16 ***
item 0.031751 0.003048 10.42 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2001 on 4451 degrees of freedom
Multiple R-squared: 0.0238, Adjusted R-squared: 0.02359
F-statistic: 108.5 on 1 and 4451 DF, p-value: < 2.2e-16
这是数据图以及 lm 和 gam 的预测,两者看起来完全合理:
plot(logart~item,ev,col="#20202020")
abline(lm(logart~item,ev),col="red")
points(ev$item,predict(me,ev),col="blue")
因此,整体模型看起来不错,只是平滑实际上并未去趋势,并且似乎完成了大部分工作(注意 y 轴范围):
plot.gam(me)
(结果与
bam
相同,并且明确指定 bs='tp'
没有区别,正如预期的那样。)
这可能吗?就与 lm/lmer 估计相匹配的 gam/bam (mis) 线性分量而言,不同的数据集如何导致如此不同的结果?我是否误解了单独测试非线性的要点(或过程)? 我还应该/还能怎么做?
以下答案由 Simon Wood 通过电子邮件发送:
平滑基础与线性趋势不正交,因此您 不会期望斜率估计值与 lm 的斜率估计值相匹配。重点 这种结构的真正目的是为了测试是否有任何东西 是否需要超越线性趋势,即是否需要平滑 而不是简单的线性效应。
我不能说我理解这个模型中的线性斜率到底代表什么,但听起来该模型不应该像我错误地假设的那样产生线性和去趋势非线性分量的分离。为此,伍德建议可以
克隆平滑构造函数并添加代码以正交化基础