library(data.table) library(dplyr) path_models = "C:/Users/antonio/ownCloud/Model_averaging/Example_AAPS/" ## get table of result files with model identifier ("model"), parameter name ("param"), parameter estimate ("value"), standard error ("SE"), ## covariances between parameters of interest ("covar") and AIC ("AIC") data_CI = data.frame(model=0,param=0,value=0,SE=0,covar_R0_delta=0,covar_R0_p=0,covar_delta_p=0,AIC=0) model = c(1:9) param = c("R0_pop","p_pop","delta_pop") loop <- txtProgressBar(min = 0, max = length(model)*length(param), style = 3) for ( k in c(1:length(model))){ for ( m in 1:length(param)){ result<-read.table(paste0(path_models,"md_",model[k],"/project/populationParameters.txt"),header=TRUE,sep=",",dec=".") # Output of estimated parameters var_covar = read.table(paste0(path_models,"md_",model[k],"/project/FisherInformation/covarianceEstimatesSA.txt"),header=F,sep=",",dec=".") BIC = read.table(paste0(path_models,"md_",model[k],"/project/LogLikelihood/logLikelihood.txt"), sep=",", header=T) # unlist(strsplit(BIC[33,],":")) rownames(result) = result[,1] data_CI[m+ (k-1)*length(param),"param"] = param[m] data_CI[m+ (k-1)*length(param),"value"] = result[param[m],"value"] data_CI[m+ (k-1)*length(param),"SE"] = result[param[m],"se_sa"] # flag_BIC = unlist(strsplit(as.character(BIC[BIC$V1=="Akaike Information Criteria (AIC)",]$V2),split=c(" "),fixed=F)) # len = length(flag_BIC) data_CI[m+ (k-1)*length(param),"covar_R0_delta"] = var_covar[1,3] data_CI[m+ (k-1)*length(param),"covar_R0_p"] = var_covar[1,4] data_CI[m+ (k-1)*length(param),"covar_delta_p"] = var_covar[2,3] data_CI[m+ (k-1)*length(param),"AIC"] = BIC[2,"importanceSampling"] data_CI[m+ (k-1)*length(param),"model"] = paste(model[k]) setTxtProgressBar(loop, m+ (k-1)*length(param)) } } ## Computation of the weights res = data_CI[data_CI$param=="R0_pop",] %>% mutate(ref = min(AIC), weight = exp(-(AIC-ref)/2)/sum(exp(-(AIC-ref)/2))) res$model = factor(res$model,levels=res$model,labels=c('k=4 & V[0]=10^3','k=4 & V[0]=10^4','k=4 & V[0]=10^5', 'k=6 & V[0]=10^3','k=6 & V[0]=10^4','k=6 & V[0]=10^5', 'k=8 & V[0]=10^3','k=8 & V[0]=10^4','k=8 & V[0]=10^5')) ## Histogram of the weights according to each candidate model. ## Figures shows that 6 models are consistent with the data ggplot(res,aes(x=model,y=weight,fill=model))+ geom_bar(stat ='identity')+ geom_text(aes(label=format(round(AIC,2), nsmall = 2),x=model,y=weight+0.01))+ theme_classic(14)+ scale_y_continuous(breaks=c(0,0.05,0.1,0.15,0.2,0.25))+ scale_x_discrete(labels=c(expression(k==4~'&'~V[0]==10^3),expression(k==4~'&'~V[0]==10^4),expression(k==4~'&'~V[0]==10^5), expression(k==6~'&'~V[0]==10^3),expression(k==6~'&'~V[0]==10^4),expression(k==6~'&'~V[0]==10^5), expression(k==8~'&'~V[0]==10^3),expression(k==8~'&'~V[0]==10^4),expression(k==8~'&'~V[0]==10^5)))+ labs(x="Candidate models",y="Weights")+ theme(axis.text.x = element_text(angle=45,hjust=1), legend.position = "none", panel.grid.major = element_line(colour="grey85")) ## Computation of the model-averaged confidence intervalss of R0, delta and p parameters model = c(1:9) nrep = 1e3 mix = NULL mix[["rep"]] = rep(c(1:1e3),times=length(model)) mix[["model"]] = rep(c(1:length(model)),each=1e3) library(mvtnorm) tir = NULL ## Draw of parameters for each candidate model for(k in 1:9){ dC = data_CI[data_CI$model==model[k],] tir[[k]] = as.data.table(rmvnorm(mean=dC$value, sigma=matrix(c(dC$SE[1]^2,dC$covar_R0_p[1],dC$covar_R0_delta[1], dC$covar_R0_p[1],dC$SE[2]^2,dC$covar_delta_p[1], dC$covar_R0_delta[1],dC$covar_delta_p[1],dC$SE[3]^2),ncol=3), n=1e3)) } tir_mix = rbindlist(tir) mix = as.data.table(mix) mix$R0 = tir_mix[,1] mix$p = tir_mix[,2] mix$delta = tir_mix[,3] mix = setorder(mix,rep,model) mix$type2 = rep(sample(x=c(1:length(model)), size=1000,replace=T, prob=res_MA$weight),each=length(model)) mix_R0 = mix[,c(1,2,3,6)] names(mix_R0)[3] = 'values' mix_R0$param = 'R0' mix_p = mix[,c(1,2,4,6)] names(mix_p)[3] = 'values' mix_p$param = 'p' mix_delta = mix[,c(1,2,5,6)] names(mix_delta)[3] = 'values' mix_delta$param = 'delta' mixx = rbind(mix_R0,mix_p,mix_delta) mix_MA = mixx[mixx$model==mixx$type2,];mix_MA$method = "MA" mix_MS = mixx[mixx$model==5,];mix_MS$method = "MS" # Model 5 is the model providing the lowest AIC mix_1 = mixx[mixx$model==1,];mix_1$method = "1" mix_2 = mixx[mixx$model==2,];mix_2$method = "2" mix_3 = mixx[mixx$model==3,];mix_3$method = "3" mix_4 = mixx[mixx$model==4,];mix_4$method = "4" mix_5 = mixx[mixx$model==5,];mix_5$method = "5" mix_6 = mixx[mixx$model==6,];mix_6$method = "6" mix_7 = mixx[mixx$model==7,];mix_7$method = "7" mix_8 = mixx[mixx$model==8,];mix_8$method = "8" mix_9 = mixx[mixx$model==9,];mix_9$method = "9" ## Table of summarizing median, 2.5th and 97.5th quantiles for each model, MS and MA CI = rbind(mix_MA,mix_MS,mix_1,mix_2,mix_3,mix_4,mix_5,mix_6,mix_7,mix_8,mix_9) %>% group_by(method,param) %>% summarize(y0 = quantile(values, 0.025,na.rm=T), y50 = median(values,na.rm=T), y975 = quantile(values, 0.975,na.rm=T)) ## Plot of confidence intervals CI$method = factor(CI$method,levels=c(1,2,3,4,5,6,7,8,9,"MS","MA"),labels=c(1,2,3,4,5,6,7,8,9,"MS","MA")) CI$param = c("delta","pi","R[0]") ggplot()+ geom_point(data=CI,aes(y=y50,x=factor(method)),shape=8,col="red")+ facet_wrap(param~.,scales="free",labeller=labeller(param=label_parsed))+ scale_x_discrete(labels=c("1"=expression(k==4~'&'~V[0]==10^3),"2"=expression(k==4~'&'~V[0]==10^4),"3"=expression(k==4~'&'~V[0]==10^5), "4"=expression(k==6~'&'~V[0]==10^3),"5"=expression(k==6~'&'~V[0]==10^4),"6"=expression(k==6~'&'~V[0]==10^5), "7"=expression(k==8~'&'~V[0]==10^3),"8"=expression(k==8~'&'~V[0]==10^4),"9"=expression(k==8~'&'~V[0]==10^5), "MS"="MS","MA"="MA"))+ geom_segment(data=CI,aes(y=y0,yend=y975,x=factor(method),xend=factor(method)))+ labs(x="Method",y="Confidence intervals")+ theme_classic(14)+ theme(axis.text.x = element_text(angle=45,hjust=1), legend.position = "none")