负二项式 GLMM 的置信区间:方差参数的不确定性
Confidence Intervals for negative binomial GLMM: uncertainty in the variance parameters
我想使用 lme4
包计算负二项式 GLMM 中预测的 95 个置信区间。但是,默认情况下没有计算预测标准误差的选项的原因我尝试:
#Package
library(lme4)
# Create a data set
set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
f2 = LETTERS[1:2], g=1:9, rep=1:15,
f3 = rnorm(10,25),
KEEP.OUT.ATTRS=FALSE)
mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)
# Create a negative binomial generalized mixed model
m.nb <- glmer.nb(y ~ f1 + f2 + poly(f3,2) + (1|g), data=dd, verbose=TRUE)
m.nb
#Compute the confidence interval for model predictions
mm<-model.matrix(~f1+f2+poly(f3,2),dd)
dd$y_est<-mm%*%fixef(m.nb) #predicted values
pvar <- diag(mm %*% tcrossprod(vcov(m.nb),mm))
CIupper <- dd$y_est+sqrt(pvar) * 1.96
CIlower <- dd$y_est-sqrt(pvar) * 1.96
但我不太确定这是否是将不确定性纳入方差参数的有效方法,尤其是使用 poly()
变量时。请问有什么帮助吗?
在上面的代码中计算模型预测的置信区间是可以的,但是它需要将逆link函数(exp()
)应用到dd$y_est
和pvar
+ sigma(m.nb)^2
:
mm<-model.matrix(~f1+f2+poly(f3,2),dd)
dd$y_est<-exp(mm%*%fixef(m.nb)) #predicted values
pvar <- diag(mm %*% tcrossprod(vcov(m.nb),mm))
CIupper <- dd$y_est+sqrt(pvar+sigma(m.nb)^2) * 1.96
CIlower <- dd$y_est-sqrt(pvar+sigma(m.nb)^2) * 1.96
我想使用 lme4
包计算负二项式 GLMM 中预测的 95 个置信区间。但是,默认情况下没有计算预测标准误差的选项的原因我尝试:
#Package
library(lme4)
# Create a data set
set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
f2 = LETTERS[1:2], g=1:9, rep=1:15,
f3 = rnorm(10,25),
KEEP.OUT.ATTRS=FALSE)
mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)
# Create a negative binomial generalized mixed model
m.nb <- glmer.nb(y ~ f1 + f2 + poly(f3,2) + (1|g), data=dd, verbose=TRUE)
m.nb
#Compute the confidence interval for model predictions
mm<-model.matrix(~f1+f2+poly(f3,2),dd)
dd$y_est<-mm%*%fixef(m.nb) #predicted values
pvar <- diag(mm %*% tcrossprod(vcov(m.nb),mm))
CIupper <- dd$y_est+sqrt(pvar) * 1.96
CIlower <- dd$y_est-sqrt(pvar) * 1.96
但我不太确定这是否是将不确定性纳入方差参数的有效方法,尤其是使用 poly()
变量时。请问有什么帮助吗?
在上面的代码中计算模型预测的置信区间是可以的,但是它需要将逆link函数(exp()
)应用到dd$y_est
和pvar
+ sigma(m.nb)^2
:
mm<-model.matrix(~f1+f2+poly(f3,2),dd)
dd$y_est<-exp(mm%*%fixef(m.nb)) #predicted values
pvar <- diag(mm %*% tcrossprod(vcov(m.nb),mm))
CIupper <- dd$y_est+sqrt(pvar+sigma(m.nb)^2) * 1.96
CIlower <- dd$y_est-sqrt(pvar+sigma(m.nb)^2) * 1.96