我正在尝试估计零膨胀康威麦克斯韦泊松混合模型的参数。我不明白为什么 GlmmTMP 函数为非零效应部分提供大约一半的值,并为零部分和色散部分提供很好的估计?
例如:- 截距的实际值是 2.5,我得到 1.21,性别女性的实际值是 1.2,我得到 0.548342 请在这种情况下帮帮我吗?谢谢
#--------Simulation from ZICOMP mix lambda---------
library(COMPoissonReg)
library(glmmTMB)
set.seed(123)
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_CMP <- 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_CMP)
Z <- model.matrix(~ 1, data = DF_CMP)
# design matrices for the fixed and random effects zero part
X_zi <- model.matrix(~ sex, data = DF_CMP)
betas <- c(2.5 , 1.2 , 2.3, -1.5) # fixed effects coefficients non-zero part
shape <- 2
gammas <- c(-1.5, 0.9) # fixed effects coefficients zero part
D11 <- 0.5 # variance of random intercepts non-zero part
# we simulate random effects
b <- rnorm(n, sd = sqrt(D11))
# linear predictor non-zero part
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF_CMP$id,drop = FALSE]))
# linear predictor zero part
eta_zi <- as.vector(X_zi %*% gammas)
DF_CMP$CMP_y <- rzicmp(n * K, lambda = exp(eta_y), nu = shape, p = plogis(eta_zi))
hist(DF_CMP$CMP_y)
#------ estimation -------------
CMPzicmpm0 = glmmTMB(CMP_y~ sex*time + (1|id) , zi= ~ sex, data = DF_CMP, family=compois)
summary(CMPzicmpm0)
> summary(CMPzicmpm0)
Family: compois ( log )
Formula: CMP_y ~ sex * time + (1 | id)
Zero inflation: ~sex
Data: DF_CMP
AIC BIC logLik deviance df.resid
4586.2 4623.7 -2285.1 4570.2 792
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
id (Intercept) 0.1328 0.3644
Number of obs: 800, groups: id, 100
Overdispersion parameter for compois family (): 0.557
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.217269 0.054297 22.42 < 2e-16 ***
sexfemale 0.548342 0.079830 6.87 6.47e-12 ***
time 1.151549 0.004384 262.70 < 2e-16 ***
sexfemale:time -0.735348 0.009247 -79.52 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.6291 0.1373 -11.866 < 2e-16 ***
sexfemale 0.9977 0.1729 5.771 7.89e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1