##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);
}
}
最近のコメント