使用GLMMadaptive绘制零膨胀、半连续数据？

Question

我正在尝试使用这里描述的effectPlotData：https://cran.r-project.org/web/packages/GLMMadaptive/vignettes/Methods_MixMod.html

但是，我正在尝试将其应用于一个模型(零膨胀半连续数据的两部分混合模型)，该模型包括线性和逻辑部分(栅栏对数正态)的随机/固定效应。我得到以下错误：'Error in Qs1，：number of dimensions不正确‘

我认为，这是因为有超过一组随机/固定效果的结果，但如果其他人遇到这个错误或可以提供建议，将不胜感激！我尝试了更改新数据框中的术语，并使用length.out尝试了几个不同的选项(尝试使用受试者的数量，然后是所有受试者的总观测值)，但每次都得到相同的错误。

下面的代码，将模型指定为m，将新数据帧指定为nDF：

m = mixed_model(Y~X, random = ~1|Subject,
                data = data_combined_temp_Fix_Num3,
                family = hurdle.lognormal,
                n_phis = 1, zi_fixed = ~X , zi_random = ~1|Subject,
                na.action = na.exclude)


nDF <- with(data_combined_temp_Fix_Num3,
            expand.grid(X = seq(min(X), max(X), length.out = 908),
                        Y = levels(Y)))


effectPlotData(m, nDF)

Answer 1

它似乎适用于下面的示例：

library("GLMMadaptive")

set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 5 # maximum follow-up time

# we constuct a data frame with the design: 
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
                 time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
                 sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects non-zero part
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
# design matrices for the fixed and random effects zero part
X_zi <- model.matrix(~ sex, data = DF)
Z_zi <- model.matrix(~ 1, data = DF)

betas <- c(-2.13, -0.25, 0.24, -0.05) # fixed effects coefficients non-zero part
sigma <- 0.5 # standard deviation error terms non-zero part
gammas <- c(-1.5, 0.5) # fixed effects coefficients zero part
D11 <- 0.5 # variance of random intercepts non-zero part
D22 <- 0.1 # variance of random slopes non-zero part
D33 <- 0.4 # variance of random intercepts zero part

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)), rnorm(n, sd = sqrt(D33)))
# linear predictor non-zero part
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, 1:2, drop = FALSE]))
# linear predictor zero part
eta_zi <- as.vector(X_zi %*% gammas + rowSums(Z_zi * b[DF$id, 3, drop = FALSE]))
# we simulate log-normal longitudinal data
DF$y <- exp(rnorm(n * K, mean = eta_y, sd = sigma))
# we set the zeros from the logistic regression
DF$y[as.logical(rbinom(n * K, size = 1, prob = plogis(eta_zi)))] <- 0

###############################################################################

km1 <- mixed_model(y ~ sex * time, random = ~ 1 | id, data = DF, 
                   family = hurdle.lognormal(),
                   zi_fixed = ~ sex)

km1


nDF <- with(DF, expand.grid(time = seq(min(time), max(time), length.out = 15),
                            sex = levels(sex)))

plot_data <- effectPlotData(km1, nDF)

library("lattice")
xyplot(pred + low + upp ~ time | sex, data = plot_data,
       type = "l", lty = c(1, 2, 2), col = c(2, 1, 1), lwd = 2,
       xlab = "Follow-up time", ylab = "")

local({
    km1$Funs$mu_fun <- function (eta) {
        pmax(exp(eta + 0.5 * exp(2 * km1$phis)), .Machine$double.eps)
    }
    km1$family$linkfun <- function (mu) log(mu)
    plot_data <- effectPlotData(km1, nDF)

    xyplot(exp(pred) + exp(low) + exp(upp) ~ time | sex, data = plot_data,
           type = "l", lty = c(1, 2, 2), col = c(2, 1, 1), lwd = 2,
           xlab = "Follow-up time", ylab = "")
})

Answer 2

如果有人遇到同样的错误，我在模型中过滤我的数据框中的数据--这导致模型的尺寸和数据框中的变量不匹配。我对新的数据帧应用了相同的过滤(我还向前移动了一个全新的数据帧，它只包括模型实际使用的试验，因此在任何步骤都不需要使用过滤)。

m = mixed_model(Y~X, random = ~1|Subject,
 data = data_combined_temp_Fix_Num3[data_combined_temp_Fix_Num3$Z>=4 &
 data_combined_temp_Fix_Num3$ZZ>= 4,],
 family = hurdle.lognormal,
 n_phis = 1, zi_fixed = ~X , zi_random = ~1|Subject,
 na.action = na.exclude)`


nDF <- with(data_combined_temp_Fix_Num3,
  expand.grid(X = seq(min(X[data_combined_temp_Fix_Num3$Z>= 4 &     
  data_combined_temp_Fix_Num3$ZZ>= 4])),
  max(X[data_combined_temp_Fix_Num3$Z>= 4 &
  data_combined_temp_Fix_Num3$ZZ>= 4])), length.out = 908),
  Y = levels(Y)))`


effectPlotData(m, nDF)