############################################################################### # MakeFigures.R # # Generate tables and plots for Kenya based on saved data. # # Author: Geir-Arne Fuglstad = 67)/1000 } png(paste(fRes, 'Kenya_median-MDG.png', sep = ""), height = 4.5, width = 4.5, units = "in", res = 200) plotOnProv(xMat, yMat, perCent.median, districtNum, xlab = "Longitude", ylab = "Latitude", kenyaMap = kenyaMap, cex.axis = 1.5, cex.lab = 1.5, asp = 1, col = viridis(64)) dev.off() png(paste(fRes, 'Kenya_prob-MDG.png', sep = ""), height = 4.5, width = 4.5, units = "in", res = 200) plotOnProv(xMat, yMat, prob, districtNum, xlab = "Longitude", ylab = "Latitude", kenyaMap = kenyaMap, cex.axis = 1.5, cex.lab = 1.5, asp = 1, col = viridis(64), zlim = c(0,1)) dev.off() # National level pNational = (nat1990-nat2015)/nat1990*100 print(median(pNational)) sum(pNational >= 67)/1000 ## Region estimates # Find approprate zlim medLim = c() sdLim = c() perLim = c() for(year in seq(1980, 2020, 5)){ load(paste(folderData, "Kenya_Prediction_Full_SPDE_", year, ".RData", sep = "")) medLim = c(medLim, range(u5Res$u5mr.reg.median, na.rm = TRUE)) sdLim = c(sdLim, range(u5Res$u5mr.reg.sd, na.rm = TRUE)) perLim = c(perLim, range(u5Res$u5mr.reg.sd/u5Res$u5mr.reg.med*100, na.rm = TRUE)) } medLim = range(medLim) sdLim = range(sdLim) perLim = range(perLim) perLim = c(3.8, 12) medLim[1] = 0.02 breaks = c(8, 12, 16, 24, 35, 1e100) # Plot each figure for(year in seq(1980, 2020, 5)){ load(paste(folderData, "Kenya_Prediction_Full_SPDE_", year, ".RData", sep = "")) print(range(u5Res$u5mr.reg.sd/u5Res$u5mr.reg.median)) png(paste(fRes, 'Kenya_RegionPred_year_', year, '.png', sep = ""), height = 8, width = 8, units = "in", res = 140) #pdf(paste(fRes, 'RegionPred_year_', year, '.pdf', sep = "")) ticks = c(0.02, 0.04, 0.08, 0.16, 0.27) plotHatch(regPred = list(med = u5Res$u5mr.reg.median, sd = u5Res$u5mr.reg.sd), zLim = medLim, sdLim = perLim, breaks = breaks, ticks = ticks, regMap = kenyaMap) dev.off() } ## National estimates # File name fName = "Kenya_Prediction_Full_SPDE_" # Initialize storage fullQuant = matrix(0, nrow = 41, ncol = 3) # Read from disk for(i in 1:41){ # Load full prediction load(paste(folderData, fName, 1979+i, ".RData", sep = "")) fullQuant[i,] = quantile(u5Res$national.samples$regVal[[1]], probs = c(0.025, 0.5, 0.975)) } # Read childmortality.org data UNB3 = read.csv(file = "../../Data/Kenya/Kenya_National/KenyaB3.csv", header = TRUE, sep = ";") DHS2014 = read.csv(file = "../../Data/Kenya/Kenya_National/DHS2014_direct.csv", header = TRUE, sep = ";") DHS2014mm = read.csv(file = "../../Data/Kenya/Kenya_National/DHS2014_MM_direct.csv", header = TRUE, sep = ";") DHS2009 = read.csv(file = "../../Data/Kenya/Kenya_National/DHS2009_direct.csv", header = TRUE, sep = ";") DHS2009mm = read.csv(file = "../../Data/Kenya/Kenya_National/DHS2009_MM_direct.csv", header = TRUE, sep = ";") DHS2003 = read.csv(file = "../../Data/Kenya/Kenya_National/DHS2003_direct.csv", header = TRUE, sep = ";") DHS2003mm = read.csv(file = "../../Data/Kenya/Kenya_National/DHS2003_MM_direct.csv", header = TRUE, sep = ";") # Make figure for full model compared to DHSs pdf(paste(fRes, "Kenya_NationalEstimates.pdf", sep = "")) matplot(1980:2020, fullQuant[,2], lty = 1, col = 1, type = "l", ylim = c(0.025, 0.13), lwd = 1.6, ylab = "National U5MR", xlab = "Year", xaxt = "n", cex.axis = 1.5, cex.lab = 1.5) axis(1, at=seq(1980, 2020, by=5), labels=seq(1980, 2020, by=5), cex.axis = 1.5, cex.lab = 1.5) polygon(c(1980:2020, 2020:1980), c(fullQuant[,1],rev(fullQuant[,3])), col=t_col("black"), border=NA) lines(UNB3$Year, UNB3$median/1000, lty = 1, lwd = 1.5, col = "red") polygon(c(UNB3$Year, rev(UNB3$Year)), c(UNB3$lower/1000,rev(UNB3$upper/1000)), col=t_col("red"), border=NA) lines(DHS2003$Year, DHS2003$U5MR/1000, type = "b", col = "blue", pch = 16) lines(DHS2003mm$Year, DHS2003mm$U5MR/1000, type = "b", col = "blue", pch = 17) lines(DHS2009$Year, DHS2009$U5MR/1000, type = "b", col = "green", pch = 16) lines(DHS2009mm$Year, DHS2009mm$U5MR/1000, type = "b", col = "green", pch = 17) lines(DHS2014$Year, DHS2014$U5MR/1000, type = "b", col = "purple", pch = 16) lines(DHS2014mm$Year, DHS2014mm$U5MR/1000, type = "b", col = "purple", pch = 17) legend('topright', legend = c("Model", "B3", "DHS2003", "DHS2003mm", "DHS2009", "DHS2009mm", "DHS2014", "DHS2014mm"), lty = c(1,1,1,1,1,1,1,1), col = c("black", "red", "blue", "blue", "green", "green", "purple", "purple"), pch = c(NA, NA, 16, 17, 16, 17, 16, 17)) dev.off() ## Compare model estimates to direct estimates # Load model estimates for 5-year periods fName = "Kenya_Period_Full_SPDE.RData" load(paste(folderData, fName, sep = "")) res.full = res.model # Load direct estimates load('../../Data/Kenya/kenya-resultsDirectEstimates-countylevelHIV.Rdata') result.direct = combinedRes perm = order(result.direct$region, result.direct$year) result.direct = result.direct[perm,c("region", "year", "var", "mean", "lower", "upper")] # Combine Surveys res.direct.full = list(med = matrix(NA, nrow = 47, ncol = 7), sd = matrix(NA, nrow = 47, ncol = 7)) for(rIdx in 1:47){ for(yIdx in 1:7){ idx = which((as.numeric(as.factor(result.direct$year)) == yIdx) & (as.numeric(as.factor(result.direct$region)) == rIdx)) if(length(idx) == 0) next res.direct.full$med[rIdx, yIdx] = as.numeric(as.character(result.direct$mean[idx])) res.direct.full$sd[rIdx, yIdx] = sqrt(as.numeric(as.character(result.direct$var[idx]))) } } # Plot against each other pdf(paste(fRes, 'Kenya_ModelVSdirect_Full_suppMat.pdf', sep = "")) plot(10, xlab = "Direct estimate (logit)", ylab = "Model estimate (logit)", xlim = c(-4.5, -0.5), ylim = c(-4.5,-0.5), cex.axis = 1.5, cex.lab = 1.5) for(i in 1:7){ points(res.direct.full$med[,i], res.full$logit.m[,i], xlim = c(-5, 0), ylim = c(-5, 0), col = i, pch = 16, cex = .7) points(res.direct.full$med[,i], res.full$logit.m[,i], xlim = c(-5, 0), ylim = c(-5, 0), col = 1, pch = 1, cex = .7) } abline(a = 0, b = 1, col = "red", lwd = 1.5) legend('topleft', legend = c("80-84", "85-89", "90-94", "95-99", "00-04", "05-09", "10-14"), col = 1:7, pch = rep(16,7), bty = "o") dev.off() pdf(paste(fRes, 'Kenya_ModelVSdirect_Full_paper.pdf', sep = "")) plot(10, xlab = "Direct estimate (logit)", ylab = "Model estimate (logit)", xlim = c(-4.5, -0.5), ylim = c(-4.5,-0.5), cex.axis = 1.5, cex.lab = 1.5) for(i in 1:7){ points(res.direct.full$med[,i], res.full$logit.m[,i], xlim = c(-5, 0), ylim = c(-5, 0), col = 1, pch = 16, cex = .7) } abline(a = 0, b = 1, col = "black", lwd = 1.5) dev.off() for(i in 1:7){ pdf(paste(fRes, 'Kenya_Model_VS_Direct_Full_Period', i, '_.pdf', sep = "")) plotCompEstimates(predA = list(mu = res.direct.full$med[,i], sd = res.direct.full$sd[,i]), predB = list(mu = res.full$logit.m[,i], sd = res.full$logit.sd[,i]), xlim = c(-5, 0), cex.lab = 1.5, cex.axis = 1.5) dev.off() } ## MSE comparison # Load hold-out results cont fName = "Kenya_Period_Hold1_SPDE.RData" load(paste(folderData, fName, sep = "")) res.hold.cont = res.model # Load hold-out result Besag fName = "Kenya_period_Hold1_Besag.RData" load(paste(folderData, fName, sep = "")) res.hold.besag = res.model # Load direct estimates load('../../Data/Kenya/kenya-resultsDirectEstimates-countylevelHIVfirstFold2014.Rdata') result.direct.hold = combinedRes perm = order(result.direct.hold$region, result.direct.hold$year) result.direct.hold = result.direct.hold[perm,c("region", "year", "var", "lower", "upper", "mean")] # Load survey into matrix res.direct.hold = list(med = matrix(NA, nrow = 47, ncol = 7), sd = matrix(NA, nrow = 47, ncol = 7)) for(rIdx in 1:47){ for(yIdx in 1:7){ idx = which((as.numeric(as.factor(result.direct.hold$year)) == yIdx) & (as.numeric(as.factor(result.direct.hold$region)) == rIdx)) if(length(idx) == 0) next res.direct.hold$med[rIdx, yIdx] = as.numeric(as.character(result.direct.hold$mean[idx])) res.direct.hold$sd[rIdx, yIdx] = sqrt(as.numeric(as.character(result.direct.hold$var[idx]))) } } # Load "truth for comparison load('../../Data/Kenya/kenya-resultsDirectEstimates-countylevelHIVholdOutFirstFold.Rdata') result.direct.true = dat perm = order(result.direct.true$region, result.direct.true$year) result.direct.true = result.direct.true[perm,c("survey", "region", "years", "var.est", "lower", "upper", "logit.est.hiv")] # Load survey into matrix res.direct.true = list(med = matrix(NA, nrow = 47, ncol = 7), sd = matrix(NA, nrow = 47, ncol = 7)) for(rIdx in 1:47){ for(yIdx in 1:7){ idx = which((as.numeric(as.factor(result.direct.true$year)) == yIdx) & (as.numeric(as.factor(result.direct.true$region)) == rIdx) & (result.direct.true$survey == 2014)) if(length(idx) == 0) next res.direct.true$med[rIdx, yIdx] = as.numeric(as.character(result.direct.true$logit.est.hiv[idx])) res.direct.true$sd[rIdx, yIdx] = sqrt(as.numeric(as.character(result.direct.true$var.est[idx]))) } } # Helper function expit = function(x){ return(exp(x)/(1+exp(x))) } # Calculate mean square errors residual.model.cont = res.hold.cont$logit.m[,1:7]-res.direct.true$med residual.model.besag = res.hold.besag$logit.m[,1:7]-res.direct.true$med residual.direct = res.direct.hold$med - res.direct.true$med res.frame = data.frame(period = c("1980-1984", "1985-1989", "1990-1994", "1995-1999", "2000-2004", "2005-2009", "2010-2014"), cont.mse = apply(residual.model.cont^2, 2, mean), besag.mse = apply(residual.model.besag^2, 2, mean), direct.mse = NA, direct.miss = NA) for(i in 3:7){ idx = which(residual.direct[,i] >= -5) res.frame$direct.miss[i] = 47-length(idx) res.frame$direct.mse[i] = mean(residual.direct[idx,i]^2) } print(format(res.frame, digit = 2, scientific = TRUE)) ## Plot mesh mesh = rEffects[[3]]$dontUse$mesh pdf(paste(fRes, 'Kenya_mesh.pdf', sep = "")) plot(1000, xlim = c(30.5, 45), ylim = c(-7.5, 8.5), cex.axis = 1.5, cex.lab = 1.5, xlab = "Longitude", ylab = "Latitude", asp = 1) plot(mesh, add = TRUE) plot(kenyaMap, add = TRUE, lwd = 3) dev.off() pdf(paste(fRes, 'Kenya_grid.pdf', sep = "")) plot(1000, xlim = range(xs), ylim = range(ys), cex.axis = 1.5, cex.lab = 1.5, xlab = "Longitude", ylab = "Latitude", asp = 1) grid(nx = 300, ny = 300, lty = 1, col = "black", lwd = 0.2) plot(kenyaMap, add = TRUE, lwd = 1.5) dev.off() pdf(paste(fRes, 'Kenya_grid_zoom.pdf', sep = "")) plot(1000, xlim = c(33.9, 34.9), ylim = c(0.1,1.2), cex.axis = 1.5, cex.lab = 1.5, xlab = "Longitude", ylab = "Latitude", asp = 1) grid(nx = 30, ny = 30, lty = 1, col = "black", lwd = 0.2) plot(kenyaMap, add = TRUE, lwd = 1.5) dev.off()