########################################################
## A function that calculates a missing-data adjusted ##
## AIC due to Hens et al. (2006)                      ##
########################################################

AICW <- function(somematrix,response,covariates,model=c("lm","logistic"),nonparam=TRUE){

library(mgcv)

## preparation
somematrix <- as.data.frame(somematrix)
n <- nrow(somematrix)
m <- ncol(somematrix)
colnames(somematrix)<-paste("x",1:m, sep="")      


## missing data indicator
delta <- c(rep(NA,n))
for(i in 1:n){
if(any(is.na(somematrix[i,]))==TRUE){delta[i]<-0}
else{delta[i]<-1}
}

## check variables that we use for the calculation of missingness-probability and their "properties"
usedvar_dis <- c(NA)
usedvar_con <- c(NA)      
for(j in c(1:(dim(somematrix)[2]))){if(any(is.na(somematrix[j]))==FALSE && length(diag(ftable(somematrix[,j],somematrix[,j]))) >= 5 ){usedvar_con <- c(usedvar_con,j)}}
for(j in c(1:(dim(somematrix)[2]))){if(any(is.na(somematrix[j]))==FALSE && length(diag(ftable(somematrix[,j],somematrix[,j]))) < 5 ){usedvar_dis <- c(usedvar_dis,j)}}
usedvar_con <- usedvar_con[-c(1)]
usedvar_dis <- usedvar_dis[-c(1)]

myindependent <- paste("delta", "~", sep="")
mycovariables_con1 <- paste("s(somematrix$x", c(usedvar_con), ")", collapse="+", sep="")
mycovariables_con2 <- paste("somematrix$x", c(usedvar_con), collapse="+", sep="")
mycovariables_dis <- paste("somematrix$x", c(usedvar_dis), collapse="+", sep="")  
   
mygammodel1 <-  as.formula(paste(myindependent,mycovariables_dis,"+",mycovariables_con1))
mygammodel2 <-  as.formula(paste(myindependent,mycovariables_dis,"+",mycovariables_con2))


## missingness-probability
response2 <- somematrix[,response]
if(nonparam==TRUE){mygam <- gam(mygammodel1, family=binomial)}
if(nonparam==FALSE){mygam <- gam(mygammodel2, family=binomial)}
probabilit <- predict.gam(mygam, type="response")
weight <- delta/probabilit
weight <- as.vector(weight)
weight <- as.numeric(weight)

##calculation of my model
myindependent <- paste("somematrix$x", response, "~", sep="")
mycovariables <- paste("somematrix$x", covariates, collapse="+", sep="")
myformula <-  as.formula(paste(myindependent,mycovariables))
if(model=="lm"){mymodel <- lm(myformula, weights=weight)}
if(model=="logistic"){mymodel <- glm(myformula, family=quasibinomial, weights=weight)}

## calculation of weighted AIC
res<-residuals(mymodel)
mycases <- noquote(names(res))
weight2 <- weight[as.numeric(mycases)]
mycases2 <- noquote(names((residuals(mymodel))[residuals(mymodel)!=0]))
mu <- as.matrix(cbind(rep(1,n),somematrix[,covariates]))%*%coefficients(mymodel)
mymu   <- 1/(1+exp(-mu))
for(j in 1:length(mymu)){if(is.na(mymu[j])==TRUE){mymu[j]<-0.5}}
logmymu <- log(mymu)
logeinsminusmymu <- log(1-mymu)

if(model=="lm"){AICW<- sum(weight2)*log((weight2%*%(res^2))/(sum(weight2)))+ 2*(length(covariates)+2)}
if(model=="logistic"){AICW <- -2*(sum(weight*(response2*logmymu+(1-response2)*logeinsminusmymu)))+2*(length(covariates)+1)}

AICW<-as.vector(AICW)

return(AICW)

}