#  This is code that I wrote to run NORM using R. I cannot find the
#  document that goes with it.
library(norm)
# Make appropriate changes in file and variable names.
setwd ("c:\\Users\\Dave Desktop\\Dropbox\\Webs\\StatPages\\More_Stuff\\Missing_Data")
x <- read.table("survrateMissingNA.dat", header = TRUE)
y <- as.matrix(x)  #convert table to matrix
cat("Logistic regression using 132 cases and missing data \n\n")
print(summary(glm(formula = outcome01~survrate + gsi + avoid + intrus, binomial, data = x)))  #Use original data with missing values
attach(x)


##########
## Important  The following code will run m = 5 times. The data will be concatenated into ComFile and then analyzed.


#  Data Augmentation using norm.R
m <- 5    #Number of imputations
k <- 9    #Number of variables in raw data file
l <- 5    #Number of variables actually used in regression
CombFile <- matrix(nrow = 0, ncol = k)

for (i in 1:m) {
   s <- prelim.norm(y)  #get preliminary statistics for the analysis
   thetahat <-em.norm(s)  #Get MLE for start value
   rngseed(25672)
   theta <- da.norm(s, thetahat, steps=200, showits=TRUE) # GET MLE
   getparam.norm(s, theta) #Print out those results
   impdata <-imp.norm(s, theta, y)  #Impute the data
   filename <- paste("CombFile", i, sep = "")
   CombFile <- rbind(CombFile, impdata)
   write(t(impdata), file = "impsurvrate", ncolumns = 9, sep = " ")
   z <- data.frame(impdata)
   z$outcome01 <- round(z$outcome01, digits = 0)
   summary((glm(formula = outcome01~survrate + gsi + avoid + intrus, binomial, data = z)))  #Use imputed data.
}

## Creating the final data file with imputed data  660 rows
nPerImp <- nrow(CombFile)/m
imps <- rep(1:m, each = nPerImp)
# Add a variable representing the imputation number.
data <- as.data.frame(cbind(imps, CombFile))
data$outcome01 <- round(data$outcome01, digits = 0)
# head(data)
attach(data)

## Set up variables to hold results
b <- matrix(NA,nrow = m, ncol = 2*l)
meanb <- numeric(l)
meanvar <- numeric(l)
varb <- numeric(l)
TT <- numeric(l)
sqrtt <- numeric(l)
t <- numeric(l)

## Run a logistic regression on each of the 5 imputed data sets and store the
## coefficients and theire standard errors.
for (i in 1:m) {         # Modify following line appropriately
       model <- glm(outcome01~survrate + gsi + avoid + intrus ,subset = (imps ==i), binomial, data = data)
       a <- summary(model)
       # print(a)
       n <- 2*l
       b[i,] <-  a$coefficients[1:n]
}

##  Calculate the coefficients, st. errors, and t values across 5 imputations
for (i in 1:l) {
  meanb[i] <- mean(b[,i])
  meanvar[i] <- mean((b[,i+l]^2))
  varb[i] <- var(b[,i])
  }
  
  cat("\n\n\nThe mean regression coefficients are: \n\n")
  print(meanb)

  
for (i in 1:l) {
  TT[i] <- meanvar[i] + (1 + 1/5)*varb[i]
  sqrtt[i] <- sqrt(TT[i])
  t[i] <- meanb[i]/sqrtt[i]
  }
cat("The standard errors are: \n\n")
print(sqrtt)
cat("\n The t values are: \n\n")
print(t)