R mgcv封装变量分组实现中的广义加法混合模型GAMM

如何解决R mgcv封装变量分组实现中的广义加法混合模型GAMM

我正在尝试使用mgcv来建模非线性多级模型，其中Time为主要自变量+ 2级协变量（Xc）。由于这些数据是个人内部的，因此我希望模型能够反映1级（Person_ID）的个人内部响应。问题在于，由于我不知道该组的编码方式，我无法理解mgcv软件包的文档（请参见我的代码中的random=~1|Person_ID，这是错误的）。

以下是详细信息，数据和代码如下：因变量= DV，时间变量= Time，级别2协变量= Xc，组= Person_ID

示例数据和代码：

#Import data.
df2 <- structure(list(Person_ID = c(100003L,100003L,100006L,100016L,100027L,100031L,100032L,100033L,100052L,100060L,100078L,100079L,100088L,100099L,100100L,100106L,100122L,100129L,100138L,100140L,100142L,100147L,100159L,100166L,100167L,100185L,100187L,100197L,100207L,100208L,100247L,100250L,100270L,100271L,100272L,100273L,1000034L,1000103L,1000136L,1000143L,1000166L,1000196L,1000215L,1000226L,1000255L,1000279L,1000281L,1000282L,1000320L,1000329L,1000383L,1000401L,1000411L,1000488L,1000496L,1000535L,1000632L,1000685L,1000732L,1000735L,1000736L,1000741L,1000758L,1000821L,1000825L,1000838L,1000870L,1000880L,1000882L,1000945L,1000993L,1001010L,1001036L,1001036L
),Time = c(5L,6L,7L,2L,3L,4L,5L,1L,8L,9L,10L,11L,12L,13L,14L,5L),Xc = c(4.33333333333333,4.33333333333333,5.16666666666667,4.66666666666667,3.5,4.5,5,5.33333333333333,6,2.5,2.83333333333333,4,6.66666666666667,7,3.83333333333333,4.16666666666667,3.66666666666667,5.5,1.33333333333333,5.83333333333333,3.33333333333333,5.66666666666667,3.16666666666667,6.83333333333333,6.33333333333333,4.83333333333333,3.5),DV = c(4,2,2.25,3,NA,4.25,3.25,1.5,1,4.75,3.75,2.75,1.75,6.75,5.25,5.75,6.25,NA)),row.names = c(3L,34L,91L,90L,89L,88L,107L,106L,108L,112L,111L,110L,118L,117L,116L,119L,120L,121L,176L,177L,178L,179L,180L,181L,203L,204L,205L,257L,258L,259L,260L,261L,262L,263L,268L,269L,270L,271L,272L,273L,285L,286L,287L,288L,300L,301L,302L,303L,307L,306L,308L,309L,310L,311L,312L,313L,324L,323L,322L,333L,335L,336L,339L,338L,340L,341L,342L,343L,344L,345L,346L,351L,352L,353L,370L,369L,375L,374L,373L,372L,371L,379L,378L,377L,376L,382L,381L,383L,384L,385L,386L,387L,388L,389L,390L,391L,396L,412L,411L,410L,413L,414L,419L,418L,417L,416L,415L,439L,448L,447L,446L,458L,457L,459L,460L,461L,462L,463L,464L,465L,466L,467L,468L,469L,470L,471L,472L,473L,475L,474L,476L,477L,24L,23L,27L,26L,33L,36L,35L,47L,46L,62L,61L,64L,63L,87L,86L,95L,94L,97L,96L,98L,113L,115L,114L,141L,140L,143L,145L,144L,172L,171L,173L,187L,186L,217L,216L,228L,227L,250L,252L,251L,254L,253L,255L,256L,279L,280L,282L,281L,284L,283L,289L,291L,290L,297L,305L,304L,314L,319L,318L),class = "data.frame")
#Attempt to model.
Model1 <- gamm(DV~s(Time)+s(Xc),data=df2,random=~1|Person_ID) #Note the attempt to group by Person_ID
summary(Model1)
plot(Model1,shift=coef(Model1)[1],pages=4,all.terms=T,rug=F,residuals=F,se=T,shade=T,seWithMean=T)

解决方法

如random中gamm参数的文档中所述：

The (optional) random effects structure as specified in a call to lme: only the list form is allowed

您必须将Person_ID变成一个因素：

df2$Person_ID <- as.factor(df2$Person_ID)

，然后在参数调用中使用一个列表：

Model1 <- gamm(DV~s(Time)+s(Xc),data=df2,random=list(Person_ID=~1)) 

您还可以使用gam来为简单的随机效果建模，如下所示：

mod <- gam(DV ~ s(Time) + s(Xc) + s(Person_ID,bs = "re"),data=df2)

这将假设Person_ID具有正常先验且具有共同方差。