Calcolare gli intervalli di confidenza per i dati medi del modello usando il restringimento in R

Question

Sto provando a eseguire un modello di sopravvivenza del nido usando il metodo di esposizione logistica basato su Shaffer, 2004. Ho una serie di parametri e desidero confrontare tutti i possibili modelli e quindi stimare i parametri mediati dal modello usando il restringimento come in Burnham e Anderson, 2002. Tuttavia, sto avendo problemi a capire come stimare gli intervalli di confidenza per i parametri aggiustati per il restringimento.

È possibile stimare intervalli di confidenza per i parametri mediati dal modello stimati utilizzando il restringimento? Posso facilmente estrarre le stime medie per i parametri mediati dal modello con il restringimento utilizzando model.average $ coef.shrinkage ma non sono chiaro come ottenere gli intervalli di confidenza corrispondenti.

Qualsiasi aiuto è gradito. Attualmente sto lavorando con il pacchetto MuMIn mentre ottengo errori con AICcmodavg riguardo alla funzione di collegamento.

Qui di seguito è una versione semplificata del codice che sto utilizzando:

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) 

Tutte le idee su come estrarre/generare i corrispondenti intervalli di confidenza?

Grazie Amy

Answer 1

Dopo riflettere su questo per un po ', hanno escogitato la seguente soluzione basata su equazione 5 in Lukacs, P. M., Burnham, K. P., & Anderson, D. R. (2009). Bias di selezione del modello e paradosso di Freedman. Annali dell'Istituto di matematica statistica, 62 (1), 117-125. Sarebbe gradito qualsiasi commento sulla sua validità.

Il codice fa seguito a quello di cui sopra:

# MuMIn generated shrinkage estimate 
    shrinkage.coef <- mod.avg$coef.shrinkage 

# beta parameters for each variable/model combination 
coef.array <- mod.avg$coefArray 
    coef.array <- replace(coef.array, is.na(coef.array), 0) # replace NAs with zeros 

# generate empty dataframe for estimates 
shrinkage.estimates <- data.frame(shrinkage.coef,variance=NA) 

# calculate shrinkage-adjusted variance based on Lukacs et al, 2009 
for(i in 1:dim(coef.array)[3]){ 
    input <- data.frame(coef.array[,,i],weight=model.set$weight) 

    variance <- rep(NA,dim(input)[2]) 
    for (j in 1:dim(input)[2]){ 
    variance[j] <- input$weight[j] * (input$Std..Err[j]^2 + (input$Estimate[j] - shrinkage.estimates$shrinkage.coef[i])^2) 
    } 
    shrinkage.estimates$variance[i] <- sum(variance) 
} 

# calculate confidence intervals 
shrinkage.estimates$lci <- shrinkage.estimates$shrinkage.coef - 1.96*shrinkage.estimates$variance 
shrinkage.estimates$uci <- shrinkage.estimates$shrinkage.coef + 1.96*shrinkage.estimates$variance