具有因子平滑交互作用的GAM分布式滞后模型（按变量）

Question

我正在尝试比较植物的两个亚类在过去 60 年中的气候响应（因子变量

subgroups

，有 2 个水平）。在同一地块上生长的两个亚组的响应是根据与长期生长的偏差来衡量的 (

plant_growth

)。气候数据包括平均温度 (

tmean

) 和平均降水量 (

prec

)。 我使用

mgcv

的

gam()

制定了一个分布式滞后模型来测试假设，即植物亚群之间的气候响应不同：

climate_model <- gam(plant_growth ~ te(tmean, lag, by = subgroups) +
                                    te(prec, lag, , by = subgroups) +
                                    te(tmean, prec, lag, , by = subgroups) ,
                                  data = plant_data)

plant_data

是一个列表，其中包含

tmean

、

prec

和

lag

作为单独的数字矩阵，

subgroups

作为区分子组 A 和 B 的因子变量，给出植物的

ID

的字符变量，并将测量的数字

plant_growth

作为向量。

但是，问题是因子

by

变量不能与

plant_data

中的矩阵参数一起使用。错误消息如下所示：

Error in smoothCon(split$smooth.spec[[i]], data, knots, absorb.cons, scale.penalty = scale.penalty,  : 
  factor `by' variables can not be used with matrix arguments.

我想知道是否有一种方法可以将因子变量

subgroups

包含到分布式滞后模型中，以便可以在因子的两个级别之间进行比较。

我已经尝试为两个级别的子组运行两个单独的滞后模型。这很好用。但是，我无法真正比较两个模型的预测，因为平滑的拟合和参数不同。此外，通过这种方式，两个亚组的气候响应被视为完全独立。然而事实并非如此。

我用 Treeclim 包中的生长数据重现了我的问题：

library("treeclim") #Data library
data("muc_spruce") #Plant growth
data("muc_clim") #Climate data

#Format climate to wide
clim <- pivot_wider(muc_clim, names_from = month, values_from = c(temp,prec))

#Format the growth data and add three new groth time series
growth <- muc_spruce %>%
            select(-samp.depth) %>%
            mutate(year = as.numeric(row.names(muc_spruce))) %>%
            mutate(ID = 1) %>%
            rename("plant_growth" = "mucstd")

additional_growth <- data.frame()
for (i in c(1:3)){
  
A <- growth %>%
  mutate(plant_growth = plant_growth + runif(nrow(muc_spruce), min = 0, max = 0.5)) %>%
  mutate(ID = ID + i)

additional_growth <- rbind(additional_growth, A)
}

growth <- rbind(growth, additional_growth)

#Bring growth and climate data together
plant_data <- na.omit(left_join(growth, clim))
rm(A, growth, clim, muc_clim, muc_spruce, additional_growth, i) #clean

#Add the subgroups label
plant_data$subgroups <- as.factor(c(rep("A", nrow(plant_data)/2), rep("B", nrow(plant_data)/2)))

#Format for gam input
plant_data <- list(lag = matrix(1:12,nrow(plant_data),12,byrow=TRUE),
                       year = plant_data$year,
                       ID = plant_data$ID,
                       plant_growth = plant_data$plant_growth,
                       subgroups = as.factor(plant_data$subgroups),
                       tmean = data.matrix(plant_data[,c(4:15)]),
                       prec = data.matrix(plant_data[,c(16:27)]))

Answer 1

来自

?mgcv::linear.functional.terms

：

该机制可用于采用因子参数的随机效应平滑，通过使用技巧来创建因子的二维数组。只需创建一个因子向量，其中包含首尾相连的因子矩阵的列（列主序）。然后重置该向量的维度以创建适当的二维数组：第一个维度应该是响应数据的数量，第二个维度应该是所需因子矩阵的列数。您不能使用matrix或data.matrix来设置所需的因子水平矩阵。请参阅下面的示例：

## set up a `factor matrix'...
fac <- factor(sample(letters,n*2,replace=TRUE))
dim(fac) <- c(n,2)

您无法创建一个艰难的因子矩阵，但可以创建一个因子并在战后修改暗淡。

Answer 2

这些实现起来具有挑战性，不幸的是，

mgcv

文档并不能让学习如何实现它们变得容易。正如上面的回复所暗示的，如果您查看线性函数（

?mgcv::linear.functional.terms

和

?mgcv::gam.models

）的帮助，您将看到

by

变量被用作加权矩阵，并且得到的平滑被构造为

rowSums(L*F) 

，其中

F

是分布式滞后平滑，

L

是加权矩阵。我在最近的几个项目上花了很多时间，最后意识到我需要在公式中为每个组提供组级分布式滞后项（这里有相当详细的描述：https://） github.com/nicholasjclark/portal_VAR）。简而言之，您可以为分组因子的每个级别创建权重矩阵，以设置组级别的分布式滞后项：

weights_l1 <- weights_l2 <- matrix(1, ncol = ncol(plant_data$lag), 
                                   nrow = nrow(plant_data$lag))

对于分组因子中的每个级别，当数据中的相应观测值与该系列匹配时，其权重矩阵中的行需要有

1

，否则为 

0

weights_l1[!(plant_data$subgroups == levels(plant_data$subgroups)[1], ] <- 0
plant_data$weights_l1 <- weights_l1

weights_l2[!(plant_data$subgroups == levels(plant_data$subgroups)[2], ] <- 0
plant_data$weights_l2 <- weights_l2

现在您已准备好根据分组因子的每个级别来拟合具有不同分布滞后项的模型

climate_model <- gam(plant_growth ~ 
                       te(tmean, lag, by = weights_l1) +
                       te(tmean, lag, by = weights_l2) +
                       te(prec, lag, by = weights_l1) +
                       te(prec, lag, by = weights_l2),
                     data = plant_data)