9

2004 年の Shaffer に基づくロジスティック暴露法を使用して、巣の生存モデルを実行しようとしています。さまざまなパラメーターがあり、考えられるすべてのモデルを比較してから、Burnham と Anderson のように収縮を使用してモデル平均パラメーターを推定したいと考えています。 2002. しかし、収縮調整済みパラメータの信頼区間を推定する方法がわかりません。

収縮を使用して推定されたモデル平均パラメータの信頼区間を推定することは可能ですか? model.average$coef.shrinkage を使用して、収縮を伴うモデル平均パラメーターの平均推定値を簡単に抽出できますが、対応する信頼区間を取得する方法が不明です。

どんな助けでも感謝します。リンク機能に関してAICcmodavgでエラーが発生するため、現在MuMInパッケージで作業しています。

以下は、私が使用しているコードの簡略版です。

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)

対応する信頼区間を抽出/生成する方法に関するアイデアはありますか?

ありがとうエイミー

4

1 に答える 1

2

これについてしばらく考えた後、私はLukacs、PM、Burnham、KP、およびAnderson、DR（2009）の方程式5に基づいて次の解決策を思いつきました。モデル選択バイアスとフリードマンのパラドックス。統計数学研究所の年報、62（1）、117–125。その有効性についてのコメントをいただければ幸いです。

上記のコードは次のとおりです。

# 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
于 2013-02-26T22:45:08.680 に答える