Berechnen Sie Konfidenzintervalle für modelldurchschnittliche Daten mit Schrumpfung in R

Question

Ich versuche, ein Nistüberlebensmodell mit der logistisch-exponentialen Methode basierend auf Shaffer, 2004, auszuführen. Ich habe eine Reihe von Parametern und möchte alle möglichen Modelle vergleichen und dann die modellgewichteten Parameter unter Verwendung von Schrumpfung schätzen, wie in Burnham und Anderson, 2002. Allerdings habe ich Schwierigkeiten herauszufinden, wie man die Konfidenzintervalle für die schrumpfungsangepassten Parameter schätzen kann.

Ist es möglich, Konfidenzintervalle für die mit Schrumpfung geschätzten modellgewichteten Parameter zu schätzen? Ich kann leicht die Mittelschätzungen für die mit Schrumpfung modellgewichteten Parameter mit Hilfe von model.average $ coef.shrinkage extrahieren, aber ich bin mir nicht sicher, wie ich die entsprechenden Konfidenzintervalle erhalten kann.

Jede Hilfe wird dankbar geschätzt. Ich arbeite derzeit mit dem MuMIn-Paket, da ich bei AICcmodavg Fehler hinsichtlich der Link-Funktion bekomme.

Im Folgenden eine vereinfachte Version des Codes, den ich verwende:

library(MuMIn)

# Logistische Expoisitions-Link-Funktion
# Siehe Shaffer, T. 2004. Ein vereinheitlichender Ansatz zur Analyse des Nestüberlebens.
# 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")
}

# Daten zufällig generieren
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

# Globales logistisch-exponentielles Modell ausführen
glm.logexp <- glm(chick/egg ~ density * height, family=binomial(logexp(days=nest.data$exposure)), data=nest.data)

# Alle möglichen Modelle bewerten
model.set <- dredge(glm.logexp)

# Modell-gewichteter 95%-Konfidenzbereich und Schätzung der Parameter unter Verwendung von Schrumpfung
mod.avg <- model.avg(model.set, beta=TRUE)
(mod.avg$coef.shrinkage)

Ideen, wie man die entsprechenden Konfidenzintervalle extrahieren/generieren kann?

Danke Amy

Answer 1

Nachdem ich eine Weile darüber nachgedacht habe, bin ich zu folgender Lösung gekommen, die auf Gleichung 5 in Lukacs, P. M., Burnham, K. P. & Anderson, D. R. (2009) basiert. Auswahlmodell-Bias und Freedmans Paradoxon. Annals of the Institute of Statistical Mathematics, 62(1), 117–125. Kommentare zur Gültigkeit sind willkommen.

Der Code folgt dem obenstehenden:

# Von MuMIn generierte Schätzverkleinerung
shrinkage.coef <- mod.avg$coef.shrinkage 

# Beta-Parameter für jede Variable/Modellkombination
coef.array <- mod.avg$coefArray
  coef.array <- ersetzen(coef.array, is.na(coef.array), 0) # Ersetze NAs durch Nullen

# Generiere leeres Datenframe für Schätzungen
shrinkage.estimates <- Datenrahmen(shrinkage.coef, varianz=NA)

# Berechne schrumpfungsangepasste Varianz basierend auf Lukacs u.a., 2009
for(i in 1:dim(coef.array)[3]){
  eingabe <- Datenrahmen(coef.array[,,i], gewicht=model.set$gewicht)

  varianz <- rep(NA, dim(eingabe)[2])
  for (j in 1:dim(eingabe)[2]){
    varianz[j] <- eingabe$gewicht[j] * (eingabe$Std..Err[j]^2 + (eingabe$Schätzung[j] - shrinkage.estimates$shrinkage.coef[i])^2)
  }
  shrinkage.estimates$varianz[i] <- sum(varianz)  
}

# Berechne Konfidenzintervalle
shrinkage.estimates$lci <- shrinkage.estimates$shrinkage.coef - 1.96*shrinkage.estimates$varianz
shrinkage.estimates$uci <- shrinkage.estimates$shrinkage.coef + 1.96*shrinkage.estimates$varianz