#################################################
## The function for a imputation  based        ##
## on a generalzation of regression imputation ##
#################################################

imputation <- function(dataframe,method=c("gam", "meanimp", "medianimp", "general.regression"),binom.data=NULL,count.data=NULL,multinom.data=NULL, spn=TRUE){        

ptm <- proc.time()          # for calculation of imputation time
method <- match.arg(method) # method of imputation

sv="min"
library(mgcv)   # for the GAM models 
library(nnet)   # for multinomial logit models                                   

# specify dimensions and properties

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

#if we have a discrete subset, itīs time to specify it
matbinom <- binom.data
matcount <- count.data
matmultnom <- multinom.data
if(is.null(multinom.data)==FALSE){
for(r in 1:length(matmultnom)){misdata[,matmultnom[r]]<-as.factor(misdata[,matmultnom[r]])}
}
attach(misdata)   

######################
## useful functions ##
######################

## function that takes only the non-missing values of a vector
###
 completev <- function(myvector){

 counter<-0
 for(i in 1:length(myvector)){
 if(is.na(myvector[i])==TRUE){counter<-counter+1}
 }


 removemissing <- c(rep(NA,length(counter)))
 k<-1

 for(j in 1:length(myvector)){
 if(is.na(myvector[j])==TRUE){removemissing[k]<-j}
  if(is.na(myvector[j])==TRUE){k<-k+1}
 }

 myvector <- myvector[-removemissing]
 return(myvector)
 }
 
## function to check how many cases per variable are missing
###
colmis <- function(mymatrix){ 
countercolumn <- c(rep(0,dim(mymatrix)[2]))    
for(i in 1:dim(mymatrix)[1]){
for(j in 1:dim(mymatrix)[2]){
if(is.na(mymatrix[i,j])==TRUE){countercolumn[j]<-countercolumn[j]+1}
}}
return(countercolumn)
}

## now only for a vector
###
colmis2 <- function(myvector){
countercolumn <- c(rep(0,length(myvector)))
for(i in 1:length(myvector)){
if(is.na(myvector[i])==TRUE){countercolumn[i]<-1}
} 
return(countercolumn)
}

## function to check how many variables are having missing values
###
amountmissing <- function(mymatrix){
amount<-colmis(mymatrix)
k<-0
for(i in 1:length(amount)){
if(amount[i] >= 1){k<-k+1}
}
return(k)
}

## function that checks if all values in a row are missing 
###
totalrowmis <- function(mymatrix){
counterrow <- c(rep(0,dim(mymatrix)[1])) 
for(i in 1:dim(mymatrix)[2]){
for(j in 1:dim(mymatrix)[1]){
if(is.na(mymatrix[j,i])==TRUE){counterrow[j]<-counterrow[j]+1}
}}
return(counterrow)
}

## function that specifies EXACTLY the amount of missing cases per variable
###
misexact <- function(mymatrix){

counting<-function(myrow){
sum(is.na(myrow)==TRUE)
}
counterrow<-as.vector(apply(mymatrix,1,counting))

D1 <- matrix(0,nrow=max(counterrow),ncol=m)


for(i in 1:max(counterrow)){
for(j in 1:m){
for(k in 1:dim(mymatrix)[1]){
         if(counterrow[k]==i){D1[i,j]<-D1[i,j]+colmis2(mymatrix[k,])[j]}
}}}

return(D1)

}

## help function(s) to let us know with which variable we start and how we proceed
## it is especially needed to specify the formula in the gam-model
###
misexact2 <- function(mymatrix){
decr <- FALSE
if(sv=="max"){decr <- TRUE}
 
D1 <- misexact(mymatrix)
D2 <- matrix(0,nrow=nrow(D1),ncol=ncol(D1))
D3 <- matrix(0,nrow=nrow(D1),ncol=ncol(D1))

for(i in 1:dim(D2)[1]){
for(j in 1:dim(D2)[2]){

  # no binomial variables -> formula consists only of smoothterms
  if(is.null(matbinom)==TRUE){
  D2[i,j]  <-  paste("s(x",(sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j],")", sep="")}

  # binomial variables -> adjust formula
  if(is.null(matbinom)==FALSE){

  isdiscrete <- FALSE
  for(k in 1:length(matbinom)){
  if((sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j]==matbinom[k]){isdiscrete<-TRUE}}

  if(isdiscrete==FALSE){
  D2[i,j]  <-  paste("s(x",(sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j],")", sep="")}
  if(isdiscrete==TRUE){
  D2[i,j]  <-  paste("x",(sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j], sep="")}
  }
  
  # multinomial variables -> adjust formula
  if(is.null(matmultnom)==FALSE){

  isdiscrete <- FALSE
  for(k in 1:length(matmultnom)){
  if((sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j]==matmultnom[k]){isdiscrete<-TRUE}}

  if(isdiscrete==TRUE){
  D2[i,j]  <-  paste("x",(sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j], sep="")}
  }
  

}} 

return(D2)
}

## similar to `misexact2'; for y-Term
###
misexact3 <- function(mymatrix){
decr <- FALSE
if(sv=="max"){decr <- TRUE} 
D1 <- misexact(mymatrix)
D3 <- matrix(0,nrow=nrow(D1),ncol=ncol(D1))

for(i in 1:dim(D3)[1]){
for(j in 1:dim(D3)[2]){

  D3[i,j]  <-  paste("x",(sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j], sep="")

}} 
return(D3)
}

##  similar to `misexact2'; for x-Terms
###
misexact4 <- function(mymatrix){
decr <- FALSE
if(sv=="max"){decr <- TRUE} 
D1 <- misexact(mymatrix)
D4 <- matrix(0,nrow=nrow(D1),ncol=ncol(D1))

for(i in 1:dim(D4)[1]){
for(j in 1:dim(D4)[2]){

  # no discrete variables -> formula consists only of smoothterms
  D4[i,j]  <-  paste("x",(sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j],sep="")

}} 
return(D4)
}

##  similar to misexact 2; needed to specify link-function in GAM algorithm
###
misexact5 <- function(mymatrix){
decr <- FALSE
if(sv=="max"){decr <- TRUE} 
D1 <- misexact(mymatrix)
D5 <- matrix(0,nrow=nrow(D1),ncol=ncol(D1))

for(i in 1:dim(D5)[1]){
for(j in 1:dim(D5)[2]){

  # no discrete variables -> formula consists only of smoothterms
  D5[i,j]  <-  paste((sort(D1[i,], index=TRUE, decreasing=decr)$ix)[j],sep="")

}} 
return(D5)
}

## count the amount of missing data
###
countmissing <- function(mymatrix){
    counter  <- 0
    for(i in 1:dim(mymatrix)[1]){
    for(j in 1:dim(mymatrix)[2]){
    if(is.na(mymatrix[i,j])==TRUE){counter <- counter+1}
   
}}
  return(counter)
}

 

###################################################################
######################## Imputation methods #######################
###################################################################

#######################
### mean imputation ###
#######################

meanimp <- function(mydata){
for(i in 1:n){
   for(j in 1:m){
     if(is.na(mydata[i,j])==TRUE){mydata[i,j] <- mean(completev(mydata[,j]))}
  }
  }
  return(mydata)
}

#########################
### median imputation ###
#########################

medianimp <- function(mydata){
for(i in 1:n){
   for(j in 1:m){
     if(is.na(mydata[i,j])==TRUE){mydata[i,j] <- median(completev(mydata[,j]))}
  }
  }
  return(mydata)
}

######################
### GAM IMPUTATION ###
######################


gamimp <- function(mydata){

D1 <- misexact(mydata)
D2 <- misexact2(mydata)
D3 <- misexact3(mydata)
D4 <- misexact4(mydata)
D5 <- misexact5(mydata)

# check if we have rows where all values missing (then: median-imp in first entry)
totalrowmissing <- totalrowmis(misdata)
poormedian      <- median(completev(misdata[,1]))
for (z in 1:n){if(totalrowmissing[z]==m){misdata[z,1]<-poormedian}}
for (z in 1:n){if(totalrowmissing[z]==m){print("Warning: One row consists only of missing values. Only here `Median imputation' is used")}}

# start the loop(s) :-/
while(countmissing(misdata)>0){

  # consider the data due to its missing values
  for(i in 1:dim(D2)[1]){
  for(j in 1:dim(D2)[2]){
  
    # consider the data itself   
        l<-as.numeric(noquote(D5[i,j]))                                 
        if(any(is.na(misdata[,l]))==TRUE){   
        for(k in 1:dim(misdata)[1]){
         
    
        # the imputation                            
        if(is.na(misdata[k,l])==TRUE){
        if(amountmissing(misdata[k,])==i){ 
          
          # the model
          if(method=="gam"){
          mymissing     <- amountmissing(misdata[k,])
          mymissing2    <- misexact2(misdata[k,])
          mymissing3    <- misexact3(misdata[k,])
          mycovariables <- paste(mymissing2[i,c(1:(dim(D2)[2]-mymissing))], sep="")
          myindependent <- paste(D3[i,j], "~")
          myfamily <- paste("family=gaussian")
          if(is.element(D5[i,j],binom.data)==TRUE){myfamily <- paste("family=binomial")}
          if(is.element(D5[i,j],count.data)==TRUE){myfamily <- paste("family=poisson")}
          if(is.element(D5[i,j],multinom.data)==TRUE){mycovariables <- paste(mymissing3[i,c(1:(dim(D2)[2]-mymissing))], sep="")}
          myformula <-  as.formula(paste(myindependent, paste(mycovariables, collapse="+") ))
          if(is.element(D5[i,j],multinom.data)==FALSE){mygam <- gam(myformula, myfamily)}
          if(is.element(D5[i,j],multinom.data)==TRUE){mygam <- multinom(myformula, trace=FALSE)}
          newd <- misdata[k,]    
          if(is.element(D5[i,j],multinom.data)==FALSE){misdata[k,l] <- predict.gam(mygam,newdata=newd,type="response")} 
          if(is.element(D5[i,j],multinom.data)==TRUE){misdata[k,l] <- predict(mygam,newdata=newd)}
          }
          
          if(method=="general.regression"){
          mymissing     <- amountmissing(misdata[k,])
          mymissing3    <- misexact3(misdata[k,])
          mycovariables <- paste(mymissing3[i,c(1:(dim(D2)[2]-mymissing))], sep="")
          myindependent <- paste(D3[i,j], "~")
          myfamily <- noquote(paste("family=gaussian"))
          if(is.element(D5[i,j],binom.data)==TRUE){myfamily <- noquote(paste("family=binomial"))}
          if(is.element(D5[i,j],count.data)==TRUE){myfamily <- noquote(paste("family=poisson"))}
          myformula <-  as.formula(paste(myindependent, paste(mycovariables, collapse="+") ))
          if(is.element(D5[i,j],multinom.data)==FALSE){mygam <- gam(myformula, myfamily)}
          if(is.element(D5[i,j],multinom.data)==TRUE){mygam <- multinom(myformula, trace=FALSE)}
          newd <- misdata[k,]                                   
          if(is.element(D5[i,j],multinom.data)==FALSE){misdata[k,l] <- predict(mygam,newdata=newd,type="response")} 
          if(is.element(D5[i,j],multinom.data)==TRUE){misdata[k,l] <- predict(mygam,newdata=newd)}
          }
          
          
        }}        
      }} 
  }}
}

mydata<-misdata

#for count and binary-data we want only integers and no digits
#spn = Senseful Prediction of Nominal data
if(spn==TRUE){ 
if(is.null(count.data)==FALSE){mydata[,matcount]<-round(mydata[,matcount],digits=0)}
if(is.null(binom.data)==FALSE){mydata[,matbinom]<-round(mydata[,matbinom],digits=0)} }

return(mydata)

}   

####################################################################  
####################################################################

# Imputation as you have chosen
if(method=="meanimp"){misdata <- meanimp(misdata)}
if(method=="medianimp"){misdata <- medianimp(misdata)}
if(method=="gam"){misdata <- gamimp(misdata)}
if(method=="general.regression"){misdata <- gamimp(misdata)}

####################################################################
####################################################################

# final adjustation of data
if(is.null(multinom.data)==FALSE){
for(r in 1:length(matmultnom)){misdata[,matmultnom[r]]<-as.numeric(as.vector(misdata[,matmultnom[r]]))}
}

# How long did it take?
ptm2 <-  proc.time()
simulationsdauer <- ptm2-ptm
simulationsdauer <- round(simulationsdauer[1], digits=4)
finaltime <- paste("The imputation time was", simulationsdauer, "second(s)")
print(finaltime)

# This is our new Matrix
detach(misdata)
return(misdata)


#rm(list=ls(all=TRUE))
# end of function
}