R/predictive_likelihood_distr_compare.R

In davewutchiett/nets_demo_spatial_dynamics: Predicting Change and Modeling Relationships In NYC Local Retail and Population.

library(plyr)
# Kolmogorov-Smirnov tests for equality of probability distributions for predictive model likelihood on withheld empirical distribution
# model fit procedure ordering: 1: unbiased observed, 2: unbiased min dev, 3: unbiased avg, 4: penalized observed, 5: penalized min dev
kstests = function(set = 'pa'){
  if (set == 'pa'){
    boot.validation.pa = readRDS('code/pavalidation.rds')
    valids = ldply(boot.validation.pa, function(x) rbind(x$t))
  } else if (set == 'hf') {
    boot.validation.hf = readRDS('code/hfvalidation.rds')
    valids = ldply(boot.validation.hf, function(x) rbind(x$t))
  } else if (set == 'pop') {
    boot.validation.pop = readRDS('code/popvalidation.rds')
    valids = ldply(boot.validation.pop, function(x) rbind(x$t))
  }
  t12 = ks.test(valids[,48], valids[,36]) # misclassification percentage
  t13 = ks.test(valids[,42], valids[,36])
  t14 = ks.test(valids[,6], valids[,36])
  t15 = ks.test(valids[,18], valids[,36])
  t23 = ks.test(valids[,48], valids[,42])
  t24 = ks.test(valids[,48], valids[,6])
  t25 = ks.test(valids[,48], valids[,18])
  t34 = ks.test(valids[,42], valids[,6])
  t35 = ks.test(valids[,42], valids[,18])
  t45 = ks.test(valids[,6], valids[,18])
  out.table.1 = rbind(c(t12$statistic, t12$p.value, t13$statistic, t13$p.value, t14$statistic, t14$p.value, t15$statistic, t15$p.value),
    c('','', t23$statistic, t23$p.value, t24$statistic, t24$p.value, t25$statistic, t25$p.value),
    c('','', '', '', t34$statistic, t34$p.value, t35$statistic, t35$p.value),
    c('', '', '', '', '', '', t45$statistic, t45$p.value))                
  tt12 = ks.test(valids[,44], valids[,32]) #residual deviance
  tt13 = ks.test(valids[,38], valids[,32])
  tt14 = ks.test(valids[,2], valids[,32])
  tt15 = ks.test(valids[,14], valids[,32])
  tt23 = ks.test(valids[,44], valids[,38])
  tt24 = ks.test(valids[,44], valids[,2])
  tt25 = ks.test(valids[,44], valids[,14])
  tt34 = ks.test(valids[,38], valids[,2])
  tt35 = ks.test(valids[,38], valids[,14])
  tt45 = ks.test(valids[,2], valids[,14])
  out.table.2 = rbind(c(tt12$statistic, tt12$p.value, tt13$statistic, tt13$p.value, tt14$statistic, tt14$p.value, tt15$statistic, tt15$p.value),
    c('','', tt23$statistic, tt23$p.value, tt24$statistic, tt24$p.value, tt25$statistic, tt25$p.value),
    c('','', '', '', tt34$statistic, tt34$p.value, tt35$statistic, tt35$p.value),
    c('', '', '', '', '', '', tt45$statistic, tt45$p.value))
  write.csv(out.table.1, paste(set,'_misclass_KStest',Sys.Date(),'.csv',sep = ''))
  write.csv(out.table.2, paste(set,'_deviance_KStest',Sys.Date(),'.csv',sep = ''))
}
kstests(set = 'pa')
kstests(set = 'hf')
kstests(set = 'pop')