R 中具有空间自相关的 GAMM
GAMM with spatial auto-correlation in R
据我了解，GAMM存在随机误差或空间自相关误差结构。
我正在尝试 运行 具有空间自相关误差结构的 GAMM 模型，例如 corExp（参见 https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/corClasses.html）。
但我对从 GAM 建模 GAMM 感到困惑。
下面是我的 GAM（广义加性模型）代码。
library(mgcv)
mod.gam <- gam(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),data = data)
summary(mod.gam)
fits = predict(mod.gam, newdata=data, type='response')
plot(data$Chave, fits, col='blue', ylab = "Predicted AGB KNN", xlab = "Estimated AGB")
plot(mod.gam)
如何将其更改为 GAMM？当我尝试使用 GAMM（如下所述）时，我没有看到太大的变化：
model <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) +
s(BaCo_5_2) + s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),
data = data, correlation=corExp(1,form = ~ Latitude + Longitude)))
summary(model$gam)
plot(model$gam,pages=1)
据我了解，空间自相关函数应该将纬度和经度作为输入。我做得对吗？
GAMM 的输出与 GAM 的输出有何不同？
不同模型的输出：
model.gamm <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2), data = data,
correlation=corSpher(form = ~ Latitude + Longitude))
summary(model.gam$lme)
[Output]:
Linear mixed-effects model fit by maximum likelihood
Data: strip.offset(mf)
AIC BIC logLik
1176.797 1232.729 -567.3985
Correlation Structure: Spherical spatial correlation
Formula: ~Latitude + Longitude | g/g.0/g.1/g.2/g.3/g.4/g.5/g.6/g.7
Parameter estimate(s):
range
0.003845795
Fixed effects: y ~ X - 1
Value Std.Error DF t-value p-value
X(Intercept) 125.44208 5.23505 96 23.961953 0.0000
Xs(Band_3)Fx1 95.48991 38.42171 96 2.485311 0.0147
Xs(Band_7)Fx1 -17.69103 40.27462 96 -0.439260 0.6615
Xs(Band_8)Fx1 -84.45737 60.83131 96 -1.388386 0.1682
Xs(BaCo_3_2)Fx1 -183.67952 179.22609 96 -1.024848 0.3080
Xs(BaCo_5_2)Fx1 111.84971 266.21261 96 0.420152 0.6753
Xs(BaCo_5_3)Fx1 -46.55365 131.32239 96 -0.354499 0.7237
Xs(BaCo_5_4)Fx1 -10.57644 9.59420 96 -1.102378 0.2731
Xs(BaCo_8_2)Fx1 76.30654 52.36900 96 1.457094 0.1484
Xs(BaCo_8A_2)Fx1 35.37155 39.34172 96 0.899085 0.3709
#########################
model.gam <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) +
s(BaCo_5_2) + s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2), data = data)
summary(model.gamm$lme)
[Output]:
Linear mixed-effects model fit by maximum likelihood
Data: strip.offset(mf)
AIC BIC logLik
1177.106 1230.375 -568.553
Fixed effects: y ~ X - 1
Value Std.Error DF t-value p-value
X(Intercept) 123.71814 4.94795 96 25.003936 0.0000
Xs(Band_3)Fx1 102.63429 38.76940 96 2.647301 0.0095
Xs(Band_7)Fx1 -26.65351 42.00784 96 -0.634489 0.5273
Xs(Band_8)Fx1 -81.96547 61.57940 96 -1.331053 0.1863
Xs(BaCo_3_2)Fx1 -312.47086 191.52784 96 -1.631464 0.1061
Xs(BaCo_5_2)Fx1 283.77313 280.17983 96 1.012825 0.3137
Xs(BaCo_5_3)Fx1 -126.82534 139.08953 96 -0.911825 0.3641
Xs(BaCo_5_4)Fx1 -11.82504 9.69256 96 -1.220012 0.2254
Xs(BaCo_8_2)Fx1 78.95786 52.86484 96 1.493580 0.1386
Xs(BaCo_8A_2)Fx1 38.03149 40.01424 96 0.950449 0.3443
##############################
mod.gam.gam <- gam(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
+ s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),data = data)
summary(mod.gam.gam)
[Output]:
Family: gaussian
Link function: identity
Formula:
Chave ~ s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 123.718 4.155 29.78 <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(Band_3) 1.000 1.000 1.125 0.2922
s(Band_7) 5.854 6.945 2.419 0.0268 *
s(Band_8) 4.130 5.210 2.989 0.0143 *
s(BaCo_3_2) 4.724 5.694 3.128 0.0113 *
s(BaCo_5_2) 2.628 3.431 1.449 0.2275
s(BaCo_5_3) 2.614 3.262 2.713 0.0530 .
s(BaCo_5_4) 4.441 5.454 1.045 0.4176
s(BaCo_8_2) 1.000 1.000 1.144 0.2881
s(BaCo_8A_2) 3.373 4.250 1.077 0.3394
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.445 Deviance explained = 60.2%
GCV = 2578.1 Scale est. = 1829.9 n = 106
我认为你叫错了 corExp()。您使用：
corExp(1, form = ~ Latitude + Longitude)
这是将指数相关函数中相关参数的值固定为1的固定值，而不是从数据中估计，这将通过使用
来实现
corExp(form = ~ Latitude + Longitude)
如果来自 GAM (gam()) 的残差空间相关性与平滑协变量正交或非线性相关，您可能也看不到模型有太大变化：相关函数进入通过响应的协方差结构建模，而不是通过均值结构建模。
如果你想从gamm()中查看包含相关函数效果的残差，那么你需要使用：
resid(model$lme, type = "normalized")
此外，如果您想比较有结构和无结构的模型，您需要在模型的 $lme 组件上使用 anova()，因此您需要改装使用 gamm() 的 gam() 模型，只需省略 correlation 部分。那么anova(mod1$lme, mod2$lme)会给你一个指数相关函数中参数的广义似然比检验。
据我了解，GAMM存在随机误差或空间自相关误差结构。
我正在尝试 运行 具有空间自相关误差结构的 GAMM 模型，例如 corExp（参见 https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/corClasses.html）。
但我对从 GAM 建模 GAMM 感到困惑。
下面是我的 GAM（广义加性模型）代码。
library(mgcv)
mod.gam <- gam(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),data = data)
summary(mod.gam)
fits = predict(mod.gam, newdata=data, type='response')
plot(data$Chave, fits, col='blue', ylab = "Predicted AGB KNN", xlab = "Estimated AGB")
plot(mod.gam)
如何将其更改为 GAMM？当我尝试使用 GAMM（如下所述）时，我没有看到太大的变化：
model <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) +
s(BaCo_5_2) + s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),
data = data, correlation=corExp(1,form = ~ Latitude + Longitude)))
summary(model$gam)
plot(model$gam,pages=1)
据我了解，空间自相关函数应该将纬度和经度作为输入。我做得对吗？ GAMM 的输出与 GAM 的输出有何不同？
不同模型的输出：
model.gamm <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2), data = data,
correlation=corSpher(form = ~ Latitude + Longitude))
summary(model.gam$lme)
[Output]:
Linear mixed-effects model fit by maximum likelihood
Data: strip.offset(mf)
AIC BIC logLik
1176.797 1232.729 -567.3985
Correlation Structure: Spherical spatial correlation
Formula: ~Latitude + Longitude | g/g.0/g.1/g.2/g.3/g.4/g.5/g.6/g.7
Parameter estimate(s):
range
0.003845795
Fixed effects: y ~ X - 1
Value Std.Error DF t-value p-value
X(Intercept) 125.44208 5.23505 96 23.961953 0.0000
Xs(Band_3)Fx1 95.48991 38.42171 96 2.485311 0.0147
Xs(Band_7)Fx1 -17.69103 40.27462 96 -0.439260 0.6615
Xs(Band_8)Fx1 -84.45737 60.83131 96 -1.388386 0.1682
Xs(BaCo_3_2)Fx1 -183.67952 179.22609 96 -1.024848 0.3080
Xs(BaCo_5_2)Fx1 111.84971 266.21261 96 0.420152 0.6753
Xs(BaCo_5_3)Fx1 -46.55365 131.32239 96 -0.354499 0.7237
Xs(BaCo_5_4)Fx1 -10.57644 9.59420 96 -1.102378 0.2731
Xs(BaCo_8_2)Fx1 76.30654 52.36900 96 1.457094 0.1484
Xs(BaCo_8A_2)Fx1 35.37155 39.34172 96 0.899085 0.3709
#########################
model.gam <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) +
s(BaCo_5_2) + s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2), data = data)
summary(model.gamm$lme)
[Output]:
Linear mixed-effects model fit by maximum likelihood
Data: strip.offset(mf)
AIC BIC logLik
1177.106 1230.375 -568.553
Fixed effects: y ~ X - 1
Value Std.Error DF t-value p-value
X(Intercept) 123.71814 4.94795 96 25.003936 0.0000
Xs(Band_3)Fx1 102.63429 38.76940 96 2.647301 0.0095
Xs(Band_7)Fx1 -26.65351 42.00784 96 -0.634489 0.5273
Xs(Band_8)Fx1 -81.96547 61.57940 96 -1.331053 0.1863
Xs(BaCo_3_2)Fx1 -312.47086 191.52784 96 -1.631464 0.1061
Xs(BaCo_5_2)Fx1 283.77313 280.17983 96 1.012825 0.3137
Xs(BaCo_5_3)Fx1 -126.82534 139.08953 96 -0.911825 0.3641
Xs(BaCo_5_4)Fx1 -11.82504 9.69256 96 -1.220012 0.2254
Xs(BaCo_8_2)Fx1 78.95786 52.86484 96 1.493580 0.1386
Xs(BaCo_8A_2)Fx1 38.03149 40.01424 96 0.950449 0.3443
##############################
mod.gam.gam <- gam(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
+ s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),data = data)
summary(mod.gam.gam)
[Output]:
Family: gaussian
Link function: identity
Formula:
Chave ~ s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) +
s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 123.718 4.155 29.78 <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(Band_3) 1.000 1.000 1.125 0.2922
s(Band_7) 5.854 6.945 2.419 0.0268 *
s(Band_8) 4.130 5.210 2.989 0.0143 *
s(BaCo_3_2) 4.724 5.694 3.128 0.0113 *
s(BaCo_5_2) 2.628 3.431 1.449 0.2275
s(BaCo_5_3) 2.614 3.262 2.713 0.0530 .
s(BaCo_5_4) 4.441 5.454 1.045 0.4176
s(BaCo_8_2) 1.000 1.000 1.144 0.2881
s(BaCo_8A_2) 3.373 4.250 1.077 0.3394
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.445 Deviance explained = 60.2%
GCV = 2578.1 Scale est. = 1829.9 n = 106
我认为你叫错了 corExp()。您使用：
corExp(1, form = ~ Latitude + Longitude)
这是将指数相关函数中相关参数的值固定为1的固定值，而不是从数据中估计，这将通过使用
corExp(form = ~ Latitude + Longitude)
如果来自 GAM (gam()) 的残差空间相关性与平滑协变量正交或非线性相关，您可能也看不到模型有太大变化：相关函数进入通过响应的协方差结构建模，而不是通过均值结构建模。
如果你想从gamm()中查看包含相关函数效果的残差，那么你需要使用：
resid(model$lme, type = "normalized")
此外，如果您想比较有结构和无结构的模型，您需要在模型的 $lme 组件上使用 anova()，因此您需要改装使用 gamm() 的 gam() 模型，只需省略 correlation 部分。那么anova(mod1$lme, mod2$lme)会给你一个指数相关函数中参数的广义似然比检验。