使用R中的收缩率计算模型平均数据的置信区间

Question

我正在尝试使用基于Shaffer,2004的逻辑曝光方法运行嵌套生存模型.我有一系列参数,希望比较所有可能的模型,然后使用收缩率估算模型平均参数,如Burnham和Anderson,但是,我无法弄清楚如何估算收缩调整参数的置信区间.

是否可以估算使用收缩估计的模型平均参数的置信区间？我可以使用model.average $ coef.shrinkage轻松提取模型平均参数的平均估计值,但不清楚如何获得相应的置信区间.

感谢任何帮助.我正在使用MuMIn包,因为我在AICcmodavg上遇到有关链接功能的错误.

下面是我正在使用的代码的简化版本:

library(MuMIn)

# Logistical Exposure Link Function
# See Shaffer, T.  2004. A unifying approach to analyzing nest success. 
# Auk 121(2): 526-540.

logexp <- function(days = 1)
{
  require(MASS)
  linkfun <- function(mu) qlogis(mu^(1/days))
  linkinv <- function(eta) plogis(eta)^days
  mu.eta <- function(eta) days * plogis(eta)^(days-1) *
    .Call("logit_mu_eta", eta, PACKAGE = "stats")
  valideta <- function(eta) TRUE
  link <- paste("logexp(", days, ")", sep="")
  structure(list(linkfun = linkfun, linkinv = linkinv,
             mu.eta = mu.eta, valideta = valideta, name = link),
        class = "link-glm")
}

# randomly generate data
nest.data <- data.frame(egg=rep(1,100), chick=runif(100), exposure=trunc(rnorm(100,113,10)), density=rnorm(100,0,1), height=rnorm(100,0,1))
  nest.data$chick[nest.data$chick<=0.5] <- 0
  nest.data$chick[nest.data$chick!=0] <- 1

# run global logistic exposure model
glm.logexp <- glm(chick/egg ~  density * height, family=binomial(logexp(days=nest.data$exposure)), data=nest.data)

# evaluate all possible models
model.set <- dredge(glm.logexp)

# model average 95% confidence set and estimate parameters using shrinkage
mod.avg <- model.avg(model.set, beta=TRUE)
(mod.avg$coef.shrinkage)

关于如何提取/生成相应置信区间的任何想法？

谢谢艾米

Answer 1

经过一段时间的思考后，我根据 Lukacs, PM, Burnham, KP, & Anderson, DR (2009) 中的方程 5 提出了以下解决方案。模型选择偏差和 Freedman\xe2\x80\x99s 悖论。统计数学研究所年鉴，62(1), 117\xe2\x80\x93125。任何对其有效性的评论将不胜感激。

\n\n

代码继承上面的代码：

\n\n

# MuMIn generated shrinkage estimate  \n  shrinkage.coef <- mod.avg$coef.shrinkage \n\n# beta parameters for each variable/model combination\ncoef.array <- mod.avg$coefArray\n  coef.array <- replace(coef.array, is.na(coef.array), 0) # replace NAs with zeros\n\n# generate empty dataframe for estimates\nshrinkage.estimates <- data.frame(shrinkage.coef,variance=NA)\n\n# calculate shrinkage-adjusted variance based on Lukacs et al, 2009\nfor(i in 1:dim(coef.array)[3]){\n  input <- data.frame(coef.array[,,i],weight=model.set$weight)\n\n  variance <- rep(NA,dim(input)[2])\n  for (j in 1:dim(input)[2]){\n    variance[j] <- input$weight[j] * (input$Std..Err[j]^2 + (input$Estimate[j] - shrinkage.estimates$shrinkage.coef[i])^2)\n  }\n  shrinkage.estimates$variance[i] <- sum(variance)  \n}\n\n# calculate confidence intervals\nshrinkage.estimates$lci <- shrinkage.estimates$shrinkage.coef - 1.96*shrinkage.estimates$variance\nshrinkage.estimates$uci <- shrinkage.estimates$shrinkage.coef + 1.96*shrinkage.estimates$variance\n

\n