# A function that calculates AIC weights for a LIST of models and
# returns the corrresponding parameter estimate

burnham <- function(models, show.weights=TRUE){

myAIC <- c(rep(0,length(models)))
for(i in 1:length(models)){myAIC[i] <- AIC(models[[i]])}

myDAIC <- c(rep(0,length(models)))
for(i in 1:length(models)){myDAIC[i] <- AIC(models[[i]])-min(myAIC)}

expAIC <- c(rep(0,length(models)))
for(i in 1:length(models)){expAIC[i] <- exp(-0.5*myDAIC[i])}

weightAIC <- c(rep(0,length(models)))
for(i in 1:length(models)){weightAIC[i] <- (expAIC[i]/sum(expAIC))}

weightAIC <- round(weightAIC, digits=4)

mycol <- c(rep(NA,length(models)))
for(i in 1:length(mycol)){mycol[i]<-length(models[[i]]$coefficients)}
mycol2<-max(mycol)
coefficientmat <- matrix(c(rep(0,(length(models)*mycol2))), nrow=length(models), ncol=mycol2)
maxmodel <- models[(((sort(mycol, index=TRUE, decr=TRUE))$ix)[1])]
colnames(coefficientmat) <- names(coefficients(maxmodel[[1]]))

for(i in 1:length(models)){
 for(j in 1:length(models[[i]]$coefficients)){
  for(k in 1:mycol2){
 
if((names(models[[i]]$coefficients))[j] == (colnames(coefficientmat))[k]){
coefficientmat[i,k]<- (models[[i]]$coefficients)[j]
}

}}}

estimate <- weightAIC%*%coefficientmat
rownames(estimate)<- paste("averaged est.")
weightAIC <- matrix(weightAIC, ncol=length(models), nrow=1)
rownames(weightAIC)<-paste("Akaike weights")
colnames(weightAIC)<-paste("Model",1:length(models), sep=" ")

if(show.weights==TRUE){return(list(weightAIC,estimate))}
if(show.weights==FALSE){return(estimate)}

}