时间戳

首页 > 解决方案 > 使用 GLMMadaptive 绘制零膨胀、半连续数据？

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

问题描述

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

但是，我正在尝试将其应用于模型（用于零膨胀半连续数据的两部分混合模型），该模型包括线性和逻辑部分（障碍对数正态）的随机/固定效应。我收到以下错误：'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)

标签： r

解决方案

它似乎适用于以下示例：

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

推荐阅读