### This analysis of the language test data.
### $Revision: 1.1.1.1 $
### $Date: 2005/01/06 14:58:12 $

###Parameters:
asp.user <- "ralmond"
asp.dbname <- "ICT_ae"
## Set this to false if you are shipping graphs to a file, 
#onscreen=TRUE

sink(paste(asp.dbname,".out",sep=""))
postscript(paste(asp.dbname,".ps",sep=""))
onscreen=FALSE
par(ask=onscreen)


### Set things up
library("RMySQL")
asp.driver <- MySQL()
source("access.q")
source("datacube.r")
## fixes for some bugs in CODA.
source("summary.mcmc.list.r")

###Connect to database
asp.dbConnect()

###List available tables
asp.listModels()
models <- asp.listModels()

###Check number of chains and cycles
chains <- asp.getChains()
goodChains <- chains[chains>0]
window <- c(501, min(asp.getCycles()[goodChains,"chain_last_cycle"]))


### find student model and link models
## Student Models names always start with HigherEdICT
smindex <- pmatch("HigherEdICT",models)
sm.name <- models[smindex]
lms.name <- models[ - smindex]

###List available parameters
sm.pnames <- asp.listParameters(sm.name)
lms.pnames <- sapply(lms.name,asp.listParameters)

###Further break up SM parameters into iCT parameters and Skill parameters
### for seven subskills
iCTindex <- grep("ICT_Literacy",sm.pnames)
iCT.pnames <- sm.pnames[iCTindex]
skills.pnames <- sm.pnames [ - iCTindex]


###Load parameters
iCT.parameters <- asp.getParameter(sm.name,iCT.pnames,
                                    iterations=window,chain=goodChains)

## for now, we will hard code this for ICT example.  We need better
## accessor functions though.  
vnames(iCT.parameters) <- c("Proficient", "Advanced", "Basic", "Below Basic", "Minimal")

skills.parameters <- asp.getParameter(sm.name,skills.pnames,
                                    iterations=window,chain=goodChains)

## for now, we will hard code this for ICT example.  We need better
## accessor functions though.  
spn <- vnames(skills.parameters)
spn <- sub("_betaLogAlpha_0"," Intercept",spn)
spn <- sub("_betaLogAlpha_1"," log(STD)",spn)
spn <- sub("_betaLogAlpha_2"," log(Slope)",spn)
vnames(skills.parameters) <- spn


lms.parameters <- lapply(lms.name, function (model) {
  result <- asp.getParameter(model, lms.pnames[[model]],
                             iterations=window,chain=goodChains)
  lpn <- vnames(result)
  lpn <- sub("_betaLogAlpha_0"," Difficulty",lpn)
  lpn <- sub("_betaLogAlpha_1"," log(Skill Slope)",lpn)
  lpn <- sub("_betaLogAlpha_2"," log(Context Slope)",lpn)
  lpn <- sub("_betaLogAlpha_3"," Increment",lpn)
  vnames(result) <- lpn
  result
})
names(lms.parameters) <- lms.name

## Summaries
cat ("Overall Proficiency\n")
summary(iCT.parameters)
## Need to reduce because Dirichlet RVs are not of full rank.
gelman.diag(iCT.parameters[,1:4])
autocorr(iCT.parameters)

cat ("\n\n Sub-Proficiencies\n")
summary(skills.parameters)
gelman.diag(skills.parameters)
autocorr(skills.parameters)
lapply(lms.name, function (name) {
  cat ("\n\n Link Model ", name, "\n")
  lmp <- lms.parameters[[name]]
  print(summary(lmp))
  print(gelman.diag(lmp))
  print(autocorr(lmp))
})
  


## Plots
plot(iCT.parameters, ask=onscreen)
autocorr.plot(iCT.parameters, ask=onscreen)
gelman.plot(iCT.parameters[,1:4], ask=onscreen)
plot(skills.parameters, ask=onscreen)
autocorr.plot(skills.parameters, ask=onscreen)
gelman.plot(skills.parameters, ask=onscreen)
lapply(lms.parameters, function (lmp) {
  plot(lmp, ask=onscreen)
  autocorr.plot(lmp, ask=onscreen)
  gelman.plot(lmp, ask=onscreen)
})


dev.off()
sink()
###When done, disconnect
##asp.dbDisconnect()