指标函数，带有mgcv或gamm4中的按变量

问题描述

我对拟合具有多变量响应的广义加法混合模型感兴趣。一个说明问题的简单示例是一个模型，每个观察结果的形式为y_{i} = b_{i} + f(x) I_{i=2} + residual，其中i = 1,2。 b_{i}是一个拦截项，f(x)是一个平滑函数，而I_{i=2}是i=2的一个指标函数。换句话说，多元响应的第一个元素仅具有截距，第二个元素具有残差加平滑函数。

这里是显示示例的数据集。我的实际模型包括具有随机截距的聚类数据。

library(mgcv)
#> Loading required package: nlme
#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:nlme':
#> 
#>     collapse
#> The following objects are masked from 'package:stats':
#> 
#>     filter,lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect,setdiff,setequal,union
library(tidyr)

# Smooth function example,taken from mgcv::gamSim
f <- function(x) {
  0.2 * x^11 * (10 * (1 - x))^6 + 10 * 
    (10 * x)^3 * (1 - x)^10
}

dat <- tibble(
  item = list(1:2),# bivariate response
  x = runif(100) # predictor,100 observations
) %>% 
  unnest(cols = item) %>% 
  mutate(
    y = (item == 2) * f(x) + rnorm(nrow(.),sd = 1),item = factor(item)
  )

# Extract the second element of the sm list,corresponding to item = 2
sm <- smoothCon(s(x,by = item),data = dat,absorb.cons = TRUE)[[2]]
# Convert to mixed model problem,eigendecomposition
re <- smooth2random(sm,"")

# Every other row is zero,corresponding to item==1 in data. As wanted.
head(cbind(dat$item,re$Xf))
#>      [,1]       [,2]
#> [1,]    1  0.0000000
#> [2,]    2  0.3535099
#> [3,]    1  0.0000000
#> [4,]    2 -1.3431389
#> [5,]    1  0.0000000
#> [6,]    2 -0.1923331
head(cbind(dat$item,attr(re$rand[[1]],"Xr")))
#>      [,1]          [,2]        [,3]        [,4]         [,5]        [,6]
#> [1,]    1  0.0000000000  0.00000000  0.00000000  0.000000000  0.00000000
#> [2,]    2  0.0049739255  0.03091603  0.01459724 -0.006153573 -0.02637519
#> [3,]    1  0.0000000000  0.00000000  0.00000000  0.000000000  0.00000000
#> [4,]    2 -0.0708917420 -0.09857063 -0.02708073  0.143078701  0.04162007
#> [5,]    1  0.0000000000  0.00000000  0.00000000  0.000000000  0.00000000
#> [6,]    2 -0.0004066909 -0.02405613  0.01229317 -0.025445761 -0.02611414
#>            [,7]       [,8]        [,9]
#> [1,]  0.0000000  0.0000000  0.00000000
#> [2,] -0.1835807  0.1491414 -0.13907891
#> [3,]  0.0000000  0.0000000  0.00000000
#> [4,]  0.1404929 -0.3934312  0.23701909
#> [5,]  0.0000000  0.0000000  0.00000000
#> [6,] -0.1820218 -0.1334959 -0.09063485

# Create model frame
mf <- list(
  y = dat$y,item = dat$item,g = attr(re$rand[[1]],"g"),# "dummy" grouping variable
  X = re$Xf,# part of f(x) in penalty null space
  Z = attr(re$rand[[1]],"Xr") # part of f(x) in penalty range space
)

# Fit with nlme::lme
mod <- lme(y ~ 0 + item + X,random = list(g = pdIdnot(~ Z - 1)),data = mf)

# Convert to original parametrization
b.original <- re$trans.U %*% (re$trans.D * as.numeric(c(ranef(mod),fixef(mod)[["X"]])))

# Linear predictor matrix
grid <- tibble(x = seq(0,1,by = .01),item = 2)
Xp <- PredictMat(sm,data = grid)

# Function estimate,including intercept
f_est <- as.numeric(fixef(mod)[["item2"]] + Xp %*% b.original)

^{由reprex package（v0.3.0）于2020-09-22创建}

这很好用，但是需要很多编码。也可以计算置信带和p值，但需要更多的编码。

因此，我的问题是，有没有一种方法可以直接使用mgcv::gamm或gamm4::gamm4来拟合这样的模型？使用像mod <- gamm(y ~ s(x,data = dat)这样的by变量不起作用，因为对于s(x)来说item=1也不是零，这不是我想要的。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）