#install polychor library
library(polycor)

#read in transition probabilites from 3 attribute example
tprobs = read.table('transprobs3.txt', header = T)

#read the reliability function into the environment
growth_reliability = function(transprobs, numtraits){
  rel1 = matrix(NA, nrow = numtraits)
  rel2 = matrix(NA, nrow = numtraits)
  
  #loop over each trait
  for(a in 1:numtraits){
    probs00 = transprobs[,1+4*(a-1)]
    probs01 = transprobs[,2+4*(a-1)]
    probs10 = transprobs[,3+4*(a-1)]
    probs11 = transprobs[,4+4*(a-1)]
    all = as.matrix(cbind(probs00, probs01, probs10, probs11))
    
    #create 4x4 contingency table of experiment-re-experiment probabilities
    cell11 = probs00^2
    cell22 = probs01^2
    cell33 = probs10^2
    cell44 = probs11^2
    cell12 = probs00*probs01
    cell13 = probs00*probs10
    cell14 = probs00*probs11
    cell23 = probs01*probs10
    cell24 = probs01*probs11
    cell34 = probs10*probs11
    
    fourmat = matrix(NA, nrow = 4, ncol = 4)
    fourmat[1,1] = mean(cell11)
    fourmat[2,2] = mean(cell22)
    fourmat[3,3] = mean(cell33)
    fourmat[4,4] = mean(cell44)
    fourmat[1,2] = fourmat[2,1] = mean(cell12)
    fourmat[1,3] = fourmat[3,1] = mean(cell13)
    fourmat[1,4] = fourmat[4,1] = mean(cell14)
    fourmat[2,3] = fourmat[3,2] = mean(cell23)
    fourmat[2,4] = fourmat[4,2] = mean(cell24)
    fourmat[3,4] = fourmat[4,3] = mean(cell34)
    
    #compute polychoric correlation of the aggregate table
    rel1[a] = polychor(fourmat)
    
    #obtain most likely transition for each examinee
    pmax = matrix(NA, nrow = nrow(transprobs))
    for(j in 1:nrow(transprobs)){
      pmax[j] = max(all[j,])
    }
    
    #compute average of most likely class probabilities
    rel2[a] = mean(pmax)
  }
  
  results = as.matrix(cbind(rel1, rel2))
  
  #column 1 is metric 1, experiment-re-experiment reliability
  #column 2 is metric 2, average most likely transition probability
  colnames(results) = c("Exp-Reexp","Avg MLT")
  
  return(results)
  
}

#run function on transtion probs file
growth_reliability(tprobs, 3)