lmer研究（２）自動akaike weight プログラム２

lmer()のakaike weight自動計算プログラムのランダム効果複数指定可能版。

使い方は、

model<- lmer(Z~ A + B + C + (1|D/E) + (0 + C|F),・・・)
model.average.lmer2(model,random="(1|D/E) + (0 + C|F)")

とすれば、akaike weightを求めてくれます。
random=　にランダム要因部分を記入すれば計算してくれます。

lmer()はS4のため、step()が使えないので、この関数はずいぶん便利だと思います。
ちなみに、step()もlmer用に移植しようと試みていますが、いまのところまだできてません。step()およびstepAIC()のプログラムを眺めていてわかったのが、この関数はもともとlmer()の原型であるlme() (in nlme package)のオブジェクトもモデル選択できることを想定しているようだけど、実際試すとできない。今後対応するのだろうか？いまだにVRパッケージは更新しているようなので、そのうちやってくれるだろうか？いずれにせよ、lme()に対応するならlmer()に対応してくれた方が数倍便利だと思うので、maintainerのRipleyさんにがんばってほしいです。

negbin()のakaike weight自動化プログラムも同じ要領で作りたいんだけど、致命的なことにnegbin()は、update()が使えないためにできない。飯島さんのサイトにglmとglm.nbのakaike weight プログラム（下の方）を参考にして作ってみるかな。

一人でイチからプログラムを作る力はまだ自分にはないです。

「model_average_lmer2.r」をダウンロード

##Depends lme4, aod
#複数のランダム効果指定用
model.average.lmer2 <- function( object, criterion = "AIC",random, ... ) {
    # extract model variables
       terms <- object@terms;
       label <- attr(terms, "term.labels");
       kmax <- length( label );

    # define AIC function
  myAIC <- function(model, ...) {
    aic   <- AIC(logLik(model))
    if( criterion == "AIC" ) {
      aic[1];
    } else if (criterion == "AICc") {
      aic[2];
    }
  }

  # null model
  object.null <- update(object, paste(c(".~.-.",random),collapse="+") )
  model <- "(Intercept)";
  aic <- myAIC(object.null);
  coefficients <- c(summary(object.null)@coefs[1], rep(NA, kmax));
  coef.sd <- c(summary(object.null)@coefs[2], rep(NA, kmax));
  names(coef.sd) <- names(coefficients) <- c("(Intercept)", label);

  # updating models
  for ( n in 1:(2^kmax - 1) ) {
    m <- numeric();
    k <- n;
    for ( i in kmax:1 ) {
      m[i] <-  k %/% 2^(i - 1);
      k <- k %% 2^(i - 1);
    }
    m <- ifelse( m == 1, TRUE, FALSE);
    result <- update(object.null, paste(c(" .~.", label[m]), collapse = "+"));
    model[n+1] <- paste(label[m], collapse = "+");
    aic[n+1] <- myAIC(result);
    sd <- coef <- rep(NA, kmax + 1);
    coef[c(TRUE,m)] <- summary(result)@coefs[,1];
    sd[c(TRUE,m)] <- summary(result)@coefs[,2];
    coefficients <- rbind(coefficients, coef);
    coef.sd <- rbind(coef.sd, sd);
  }

  # calculating Akaike's weight w
  if( criterion == "AIC" ) {
    out <- data.frame(model, AIC=aic);
    out$dAIC <- aic - min(aic);
    out$w <- exp( -0.5 * out$dAIC ) / sum(exp( -0.5 * out$dAIC ));
    }
  else if (criterion == "AICc") {
    out <- data.frame(model, AICc=aic);
    out$dAICc <- aic - min(aic);
    out$w <- exp( -0.5 * out$dAICc ) / sum(exp( -0.5 * out$dAIC ));
    }

  # formatting coefficents and sds
  rownames(coef.sd) <- rownames(coefficients) <- 1:nrow(coefficients);

  # calculating z-statistics and p-value
  z <- coefficients / coef.sd;
  p <- pnorm(z);

  # averaging models
  ave.coef <- apply( coefficients * out$w, 2, sum, na.rm = TRUE );
  ave.var <- apply( ( coef.sd^2 + (t(t(coefficients) - ave.coef))^2 ) * out$w, 2, sum, na.rm = TRUE );
  ave.sd <- sqrt(ave.var);
  ave.z <- ave.coef / ave.sd;
  ave.p <- pnorm(ave.z);

  # importance of variables
  tmp <- ifelse(!is.na(coefficients[,-1]), matrix(out$w, nrow=nrow(out), ncol=ncol(out) - 1), 0);
  iv <- colSums(tmp);
  names(iv) <- label;

  # output
  if( criterion == "AIC" ) {
    o <- order(out$AIC);
    list(summary = out[o,], coefficients = coefficients[o,], sd = coef.sd[o,], p = p[o,],
       ave.coef = ave.coef, ave.sd = ave.sd, ave.p = ave.p, iv=iv);
    }
  else if (criterion == "AICc") {
    o <- order(out$AICc);
    list(summary = out[o,], coefficients = coefficients[o,], sd = coef.sd[o,], p = p[o,],
       ave.coef = ave.coef, ave.sd = ave.sd, ave.p = ave.p, iv=iv);
    }
}