# Please contact Ivy Liu (iliu@msor.vuw.ac.nz) if you have questions. # y.mat<-matrix(scan("Madelon500.txt",sep=","),nrow=500,ncol=501,byrow=T) y.mat[1:3,] y.cov=y.mat[,501] y.cov[y.cov==-1]<-0 y.mat=y.mat[,-501] y.mat[1:3,] y.mat[y.mat==0|y.mat==1]<-1 y.mat[y.mat==2|y.mat==3]<-2 y.mat[y.mat==4|y.mat==5]<-3 y.mat[y.mat==6|y.mat==7]<-4 y.mat[y.mat==8|y.mat==9]<-5 y.mat[1:3,] #Create vectors to store results: Names=c() #Model names Row.grs=c() #Number of row groups Column.grs=c() #Number of column groups Params=c() #Formula for number of parameters (including pi and kappa) nu.v=c() #The number of parameters in this case (including pi and kappa) AIC.v=c() # AIC values AICc.v=c() # AIC values BIC.v=c() # BIC values ICL.v=c() # ICL values #Proportional odds models # ------------------------------ # Without the covariate #-------------------------------- #1a1 covy=as.numeric(rep(y.cov,ncol(y.mat))) length(covy) length(y.mat) PO.ss.out=polr(as.factor(y.mat)~1) PO.ss.out$mu=PO.ss.out$zeta PO.ss.out$AIC=PO.ss.out$deviance+2*PO.ss.out$edf PO.ss.out$AICc=PO.ss.out$AIC+(2*(PO.ss.out$edf+1)*(PO.ss.out$edf+2))/(nrow(y.mat)*ncol(y.mat) - PO.ss.out$edf - 2) PO.ss.out$BIC=PO.ss.out$deviance+2*PO.ss.out$edf*log(nrow(y.mat)*ncol(y.mat)) Names=c(Names, "1a, PO.ss") Row.grs=c(Row.grs,1) Column.grs=c(Column.grs,1) Params=c(Params, "q-1") nu.v=c(nu.v,PO.ss.out$edf) AIC.v=c(AIC.v,PO.ss.out$AIC) AICc.v=c(AICc.v,PO.ss.out$AICc) BIC.v=c(BIC.v,PO.ss.out$BIC) ICL.v=c(ICL.v,NA) PO.ss.out #3a1 maxCG=12 POFM.sc1.out=list() for(CG in 3:maxCG){ POFM.sc1.out[[CG]]=POFM.sc1.F(invect=c(PO.ss.out$mu,rep(1,CG-1))) Names=c(Names,"3a1, POFM.sc1") Row.grs=c(Row.grs,1) Column.grs=c(Column.grs,CG) Params=c(Params, "q+2*C-3") # with no covariate nu.v=c(nu.v,POFM.sc1.out[[CG]]$info[5]) AIC.v=c(AIC.v,POFM.sc1.out[[CG]]$info[6]) AICc.v=c(AICc.v,POFM.sc1.out[[CG]]$info[7]) BIC.v=c(BIC.v,POFM.sc1.out[[CG]]$info[8]) ICL.v=c(ICL.v,POFM.sc1.out[[CG]]$info[9]) } POFM.sc1.out[[3]] data.frame("Model"=Names, Row.grs, Column.grs, "No. pars"=Params, "No. pars"=nu.v, "AIC"= AIC.v, "AICc"= AICc.v, "BIC"= BIC.v, "ICL"= ICL.v) # ------------------------------ # With the covariate #-------------------------------- #1a2 covy=as.numeric(rep(y.cov,ncol(y.mat))) length(covy) length(y.mat) PO.ss.out=polr(as.factor(y.mat)~covy) PO.ss.out$mu=PO.ss.out$zeta PO.ss.out$xcoef=PO.ss.out$coef PO.ss.out$AIC=PO.ss.out$deviance+2*PO.ss.out$edf PO.ss.out$AICc=PO.ss.out$AIC+(2*(PO.ss.out$edf+1)*(PO.ss.out$edf+2))/(nrow(y.mat)*ncol(y.mat) - PO.ss.out$edf - 2) PO.ss.out$BIC=PO.ss.out$deviance+2*PO.ss.out$edf*log(nrow(y.mat)*ncol(y.mat)) Names=c(Names, "1a, PO.ss") Row.grs=c(Row.grs,1) Column.grs=c(Column.grs,1) Params=c(Params, "q-1+1") nu.v=c(nu.v,PO.ss.out$edf) AIC.v=c(AIC.v,PO.ss.out$AIC) AICc.v=c(AICc.v,PO.ss.out$AICc) BIC.v=c(BIC.v,PO.ss.out$BIC) ICL.v=c(ICL.v,NA) PO.ss.out #3a2 maxCG=12 POFM.sc.out=list() for(CG in 3:maxCG){ POFM.sc.out[[CG]]=POFM.sc.F(invect=c(PO.ss.out$mu,PO.ss.out$xcoef,rep(1,CG-1))) Names=c(Names,"3a, POFM.sc") Row.grs=c(Row.grs,1) Column.grs=c(Column.grs,CG) #Params=c(Params, "q+2*C-3") # with no covariate Params=c(Params, "q+2*C-3+1") # with a covariate nu.v=c(nu.v,POFM.sc.out[[CG]]$info[5]) AIC.v=c(AIC.v,POFM.sc.out[[CG]]$info[6]) AICc.v=c(AICc.v,POFM.sc.out[[CG]]$info[7]) BIC.v=c(BIC.v,POFM.sc.out[[CG]]$info[8]) ICL.v=c(ICL.v,POFM.sc.out[[CG]]$info[9]) } POFM.sc.out[[3]] data.frame("Model"=Names, Row.grs, Column.grs, "No. pars"=Params, "No. pars"=nu.v, "AIC"= AIC.v, "AICc"= AICc.v, "BIC"= BIC.v, "ICL"= ICL.v) # ------------------------------ # With the covariate + interactions #-------------------------------- #1a3 covy=as.numeric(rep(y.cov,ncol(y.mat))) length(covy) length(y.mat) PO.ss.out=polr(as.factor(y.mat)~covy) PO.ss.out$mu=PO.ss.out$zeta PO.ss.out$xcoef=PO.ss.out$coef PO.ss.out$AIC=PO.ss.out$deviance+2*PO.ss.out$edf PO.ss.out$AICc=PO.ss.out$AIC+(2*(PO.ss.out$edf+1)*(PO.ss.out$edf+2))/(nrow(y.mat)*ncol(y.mat) - PO.ss.out$edf - 2) PO.ss.out$BIC=PO.ss.out$deviance+2*PO.ss.out$edf*log(nrow(y.mat)*ncol(y.mat)) Names=c(Names, "1a, PO.ss") Row.grs=c(Row.grs,1) Column.grs=c(Column.grs,1) Params=c(Params, "q-1+1") nu.v=c(nu.v,PO.ss.out$edf) AIC.v=c(AIC.v,PO.ss.out$AIC) AICc.v=c(AICc.v,PO.ss.out$AICc) BIC.v=c(BIC.v,PO.ss.out$BIC) ICL.v=c(ICL.v,NA) PO.ss.out #3a3 maxCG=12 POFM.sc3.out=list() for(CG in 4:maxCG){ POFM.sc3.out[[CG]]=POFM.sc3.F(invect=c(PO.ss.out$mu,PO.ss.out$xcoef,rep(1,CG-1),rep(1,CG-1))) Names=c(Names,"3a3, POFM.sc") Row.grs=c(Row.grs,1) Column.grs=c(Column.grs,CG) #Params=c(Params, "q+2*C-3") # with no covariate Params=c(Params, "q+2*C-3+1+C-1") # with a covariate +interactions nu.v=c(nu.v,POFM.sc3.out[[CG]]$info[5]) AIC.v=c(AIC.v,POFM.sc3.out[[CG]]$info[6]) AICc.v=c(AICc.v,POFM.sc3.out[[CG]]$info[7]) BIC.v=c(BIC.v,POFM.sc3.out[[CG]]$info[8]) ICL.v=c(ICL.v,POFM.sc3.out[[CG]]$info[9]) } POFM.sc3.out[[12]] data.frame("Model"=Names, Row.grs, Column.grs, "No. pars"=Params, "No. pars"=nu.v, "AIC"= AIC.v, "AICc"= AICc.v, "BIC"= BIC.v, "ICL"= ICL.v)