ms02 network aging simulation

Analysis using the first 15 simulated ms02 permuted networks. 

file : 20140713-compare-ms02,ginppi.R
#compare agine dynamcis between original and ms02 GINDIP.
# plot mortality rate ~time from simulated GINPPI aging.
# binomial aging model fitting

rm(list=ls())
setwd("~/github/mactower-network-failure-simulation/compareOrigMS02")
require(flexsurv)
getwd()
list.files()

source("lifespan.20140711.R")

ms02path = 'ms02.gindip.failures'
list.files(, path=ms02path)

ms02fullnames = NA;
for(i in 1:20){
  mypath = paste( ms02path, '/', i, '/popages', sep='' )
  myfiles = list.files(, path=mypath)
  myfullnames = paste( mypath, '/', myfiles, sep='')
  ms02fullnames = c(ms02fullnames, myfullnames)
}
ms02fullnames = ms02fullnames[-1]

origAgeFiles = list.files(path="ori_ginppit_ages")
origAgeFiles

mypattern = "p.0.8.lambda.0.01.popsiz"
ms02files= ms02fullnames[grep(mypattern, ms02fullnames)]
ms02output = data.frame(ms02files);
ms02output$mean = NA; ms02output$stddev = NA; ms02output$CV=NA;
ms02output$R=NA; ms02output$t0=NA; ms02output$n=NA;
ms02output$Rflex = NA; ms02output$Gflex = NA;

for( i in 1:length(ms02files)) {
  tb = read.csv(ms02files[i])
  summary(tb); names(tb) = c('rls'); tb$rls = round(tb$rls+0.5, digits=0)
  ms02output$mean[i]=mean(tb$rls)
  ms02output$stddev[i]= sd(tb$rls)
  ms02output$CV[i] = sd(tb$rls) / mean(tb$rls)
 
  fitGomp = flexsurvreg( formula=Surv(tb$rls) ~ 1, dist='gompertz')
  Rhat = fitGomp$res[2,1]; Ghat = fitGomp$res[1,1]
  ms02output$Rflex[i] = Rhat;
  ms02output$Gflex[i]=Ghat;
 
  #G = (n-1) / t0, so t0 = (n-1)/G
  nhat = 4; t0= (nhat-1)/Ghat; t0
  fitBinom = optim ( c(Rhat, t0, nhat),  llh.binomialMortality.single.run , lifespan=tb$rls,
                     method="L-BFGS-B", lower=c(1E-10,0.05,1E-1), upper=c(10,100,50) );
 
  ms02output[i, c("R", "t0", "n")] = fitBinom$par
 
  tb2 = calculate.mortality.rate(tb$rls)
  plot( tb2$s ~ tb2$t, main=ms02files[i])
  plot( tb2$mortality.rate ~ tb2$t , main=ms02files[i])
  plot( tb2$mortality.rate ~ tb2$t, log='y', main= paste( 'linear for Gompertz', main=ms02files[i]) )
  plot( tb2$mortality.rate ~ tb2$t, log='yx', main=paste( 'linear for Weibull',main=ms02files[i]) )
}#i file loop

ms02output

oripath = "ori_ginppit_ages"
#origAgeFiles = list.files(path= oripath)
origAgeFiles
origFile = origAgeFiles[grep(mypattern, origAgeFiles)]
tb.ori = read.csv(paste( oripath, '/',origFile[2],sep='')  )
tb.ori$rls = round( tb.ori[,1]+0.5, digit=0)

oriGomp = flexsurvreg( formula=Surv(tb.ori$rls) ~ 1, dist='gompertz')
Rori = oriGomp$res[2,1]; Gori = oriGomp$res[1,1]
nori = 4; t0= (nori-1)/Gori; t0
oriBinom = optim ( c(Rori, t0, nori),  llh.binomialMortality.single.run , lifespan=tb.ori$rls,
                   method="L-BFGS-B", lower=c(1E-10,0.05,1E-1), upper=c(10,100,50) );

cbind( mean(tb.ori$rls), median(ms02output$mean))
cbind( sd(tb.ori$rls), median(ms02output$stddev))
cbind( sd(tb.ori$rls)/mean(tb.ori$rls), median(ms02output$CV))
cbind( Rori, median(ms02output$Rflex))
cbind( oriBinom$par[1], median(ms02output$R))
cbind( oriBinom$par[2], median(ms02output$t0))
cbind( oriBinom$par[3], median(ms02output$n)) #ms02 does change n
cbind( Gori, median(ms02output$Gflex))

cv.ori =  sd(tb.ori$rls)/mean(tb.ori$rls)
hist(ms02output$CV, br=15)
arrows(cv.ori, 2.5, cv.ori, 1.1, lwd=2, col='red', xlab="t0" )

hist(ms02output$Rflex)
hist(ms02output$t0)
# arrows(oriBinom$par[2], 5, oriBinom$par[2], 1, lwd=2, col='red', xlab="t0" )

s.ori = calculate.mortality.rate(tb.ori$rls)
tb = read.csv(ms02files[6]); tb$rls = round(tb[,1]+0.5, digits=0)
s2 = calculate.mortality.rate(tb$rls)
ks.test( tb.ori$rls, tb$rls)
wilcox.test( tb.ori$rls, tb$rls)
var.test( tb.ori$rls, tb$rls)

plot( s.ori$s ~ s.ori$t, col='red',type='l')
lines(s2$s ~ s2$t)

plot( s.ori$mortality.rate ~ s.ori$t, type='l', col='red')
lines(s2$mortality.rate ~ s2$t, col='blue')

plot( s.ori$mortality.rate ~ s.ori$t, type='l', col='red', log='y')
lines(s2$mortality.rate ~ s2$t, col='blue')

plot( s.ori$mortality.rate ~ s.ori$t, type='l', col='red', log='xy')
lines(s2$mortality.rate ~ s2$t, col='blue')

quit('no')
##########END 20140713