Построение графиков с помощью GLMMadaptive для полунепрерывных данных с нулевым раздутием?

Question

Я пытаюсь использовать effectPlotData, как описано здесь: https://cran.r-project.org/web/packages/GLMMadaptive/vignettes/Methods_MixMod.html

Но я пытаюсь применить его к модели (двухкомпонентная смешанная модель для нуля -вдуваемые полунепрерывные данные), который включает случайные / фиксированные эффекты как для линейной части, так и для части logisti c (логарифмически-нормальный барьер). Я получаю следующую ошибку: «Ошибка в Qs [1,]: неправильное количество измерений»

Что, я думаю, связано с наличием более одного набора случайных / фиксированных результатов эффектов, но если кто-то еще пришел через эту ошибку или можете посоветовать, был бы признателен! Я попытался изменить термины в новом фрейме данных и попробовал несколько разных вариантов с 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

На случай, если кто-то обнаружит ту же ошибку, я отфильтровал данные из моего фрейма данных в модели, что привело к несовпадению размеров модели и переменной из фрейма данных. Я применил ту же фильтрацию к новому фрейму данных (я также продвинулся вперед с совершенно новым фреймом данных, который включает только испытания, которые фактически используются моделью, поэтому фильтрация не должна использоваться ни на одном этапе).

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)

Answer 2

Кажется, это работает в следующем примере:

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 = "")
})