R/test.600.VEME.MaC.incremental.R
In wdelva/MiceABC:
#setwd("/Users/wimdelva/Documents/MiceABC/R")
#source("simpact.wrapper.R")
#source("VEME.wrapper.R")
#source("dummy.wrapper.R")
#source("simpact.parallel.R")
#source("MaC.incremental.R")
#source("dummy.MaC.incremental.R")
#source("/user/data/gent/vsc400/vsc40070/phylo/scripts/VEME.wrapper.R")
#source("/user/data/gent/vsc400/vsc40070/phylo/scripts/simpact.parallel.R")
#source("/user/data/gent/vsc400/vsc40070/phylo/scripts/dummy.MaC.incremental.R")
#source("/Users/delvaw/Documents/MiceABC/R/VEME.wrapper.R")
source("/Users/delvaw/Documents/MiceABC/R/VEME.wrapper2.R")
source("/Users/delvaw/Documents/MiceABC/R/VEME.wrapper2.medians.R")
source("/Users/delvaw/Documents/MiceABC/R/mice.wrapper.R")
source("/Users/delvaw/Documents/MiceABC/R/00-Functions.R")
source("/Users/delvaw/Documents/MiceABC/R/simpact.parallel.R")
source("/Users/delvaw/Documents/MiceABC/R/mice.parallel.R")
source("/Users/delvaw/Documents/MiceABC/R/dummy.MaC.incremental.R")
source("/Users/delvaw/Documents/MiceABC/R/dummy.MaC.incremental.parallel.mice.R")
library(dplyr)
library(MASS)
library(splines)
library(boot)
#library(haven)
library(ggplot2)
library(GGally)
library(fitdistrplus)
library(lmtest)
library(mclust)
#library(depth)
library(pcaPP)
#library(devtools)
#install_github("wdelva/RSimpactHelp")
library(RSimpactCyan)
library(RSimpactHelper)
library(lmtest)
library(mice)
#library(miceadds)
library(parallel)
library(randtoolbox)
library(EasyABC)
library(dplyr)
library(tidyr)
library(nlme)
library(lme4)
library(boot)
library(data.table)
library(parallel)
install.packages("adegenet", dep = TRUE) # Still to be done on vsc
library(adegenet) # Still to be done on vsc
library(ade4) # Still to be done on vsc
source("http://adegenet.r-forge.r-project.org/files/patches/auxil.R")
#install.packages("devtools")
#library(devtools)
#install_github("emvolz/treedater") # Still to be done on vsc
library(treedater) # Still to be done on vsc
dummy.input.vector <- c(1.1, 0.25, 0, 3, 0.23, 0.23, # what if 1.1 becomes 1.4
45, 45, #45, 45, # what if 45 becomes 60
-0.5, 2.8, -0.2, -0.2, -2.5, -0.52, -0.05)# c(1000, 2, 3, 4)
x.offset <- length(dummy.input.vector)
n.experiments <- 80
dummy.master2 <- simpact.parallel(model = VEME.wrapper2,
actual.input.matrix = matrix(rep(dummy.input.vector, each = n.experiments), nrow = n.experiments),
seed_count = 0,
n_cluster = 8)
save(valid.dummy.master2, file = "/Users/delvaw/Documents/MiceABC/valid.dummy.master2.RData")
#####
# The output of the master model
#####
#inc.master.vector only exists for the validation version of the master model
# inc.master.vector <- dummy.master2[, 17]
# #save(inc.master.vector, file = "/Users/delvaw/Documents/MiceABC/inc.master.vector.RData")
# #save(inc.master.vector.inflated, file = "/Users/delvaw/Documents/MiceABC/inc.master.vector.inflated.RData")
#
# hist(inc.master.vector[!is.na(inc.master.vector)], 14) # The distribution of HIV incidence after 10 years
# mean(inc.master.vector) # The mean HIV incidence after 10 years: 0.01527118
# median(inc.master.vector) # The median HIV incidence after 10 years: 0.01532902
# quantile(inc.master.vector, c(0.025, 0.975)) # 0.002462209 0.026417416
head(dummy.master2)
inc.master.vector <- dummy.master2[, 17]
new.infect.vector <- dummy.master2[, 18]
recent.ratio.vector <- dummy.master2[, 19]
mean(new.infect.vector, na.rm = TRUE) # The mean HIV incidence after 10 years
median(new.infect.vector, na.rm = TRUE) # The median HIV incidence after 10 years
quantile(new.infect.vector, c(0.025, 0.975), na.rm = TRUE)
# mean.br.len.vector <- dummy.master2[, 19]
hist(new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf], 20)
hist(recent.ratio.vector[!is.na(recent.ratio.vector) & recent.ratio.vector < Inf], 20)
plot(new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf],
inc.master.vector[!is.na(new.infect.vector) & new.infect.vector < Inf])
plot(recent.ratio.vector[inc.master.vector > 0.0142 & inc.master.vector < 0.0162],
inc.master.vector[inc.master.vector > 0.0142 & inc.master.vector < 0.0162])
cor.test(inc.master.vector[!is.na(new.infect.vector) & new.infect.vector < Inf],
new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf])
cor.test(inc.master.vector[!is.na(new.infect.vector) & new.infect.vector < Inf],
new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf])
#plot(mean.br.len.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf],
# inc.master.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf])
#cor.test(mean.br.len.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf],
# inc.master.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf])
#plot(mean.br.len.vector,
# new.infect.vector)
dummy.master2 <- dummy.master2 %>%
as.data.frame() %>%
dplyr::filter(complete.cases(.))
dummy.targets.empirical <- l1median(dummy.master2)
# dummy.targets.empirical (based on 400 repeats):
# 5.79782475 4.16132170 0.63770526 2.69246630 1.21359926 -1.56183331 0.25472579 2.09719359 0.02891330 0.13671790
# 0.06816378 0.34554390 0.31086783 0.37050293 0.41706055 1.01682475 7.07415631
dummy.targets.empirical <- c(5.79782475,
4.16132170,
0.63770526,
2.69246630,
1.21359926,
-1.56183331,
0.25472579,
2.09719359,
0.02891330,
0.13671790,
0.06816378,
0.34554390,
0.31086783,
0.37050293,
0.41706055,
1.01682475,
7.07415631)
# round(colMeans(dummy.master2), 3) # For interest sake, what are the marginal means?
predictorMatrix <- (1 - diag(1, length(c(dummy.input.vector, dummy.targets.empirical)))) # This is the default matrix.
# # Let's now modify the first 15 rows of this matrix, corresponding to the indicators of predictor variables for the input variables. In brackets the values for the master model.
predictorMatrix[1:x.offset, ] <- 0 # First we "empty" the relevant rows, then we refill them.
# We are currently not allowing input variables to be predicted by other predictor variables. Only via output variables. We could change this at a later stage.
predictorMatrix[1, x.offset + c(10, 11, 17)] <- 1 # relative susceptibility in young women is predicted by HIV prevalence in young men and women, and recent infections (~ incidence)
predictorMatrix[2, x.offset + 3] <- 1 # agescale predicted by slope
predictorMatrix[3, x.offset + c(1, 3, 6)] <- 1 # mean of the person-specific age gap preferences is predicted by slope, intercept and AAD
predictorMatrix[4, x.offset + c(2, 4, 5)] <- 1 # sd of the person-specific age gap preferences is predicted by SD, WSD, BSD
predictorMatrix[5, x.offset + c(7, 8, 9, 13, 16)] <- 1 # man gamma a predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.men, exp(growthrate)
predictorMatrix[6, x.offset + c(7, 8, 9, 12, 16)] <- 1 # woman gamma a predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.women, exp(growthrate)
predictorMatrix[7, x.offset + c(7, 8, 9, 13, 16, 17)] <- 1 # man gamma b predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.men, exp(growthrate), and recent infections (~ incidence)
predictorMatrix[8, x.offset + c(7, 8, 9, 12, 16, 17)] <- 1 # woman gamma b predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.men, exp(growthrate), and recent infections (~ incidence)
predictorMatrix[9, x.offset + c(2, 4, 5, 7, 8, 14, 15, 16, 17)] <- 1 # formation.hazard.agegapry.gap_factor_x_exp is predicted by population growth, age gap variance, hiv prevalence, and recent infections (~ incidence)
predictorMatrix[10, x.offset + c(7, 8, 9, 12, 13, 16, 17)] <- 1 # baseline formation hazard predicted by HIV prevalence, cp, degree distrib. HIV prevalence, and recent infections (~ incidence)
predictorMatrix[11, x.offset + c(7, 8, 9, 12, 13, 16, 17)] <- 1 # numrel man penalty is predicted by degree distrib, cp, prev, popgrowth, and recent infections (~ incidence)
predictorMatrix[12, x.offset + c(7, 8, 9, 12, 13, 16, 17)] <- 1 # # numrel woman penalty is predicted by degree distrib, cp, prev, popgrowth, and recent infections (~ incidence)
predictorMatrix[13, x.offset + 16] <- 1 # conception.alpha_base is predicted by popgrowth
predictorMatrix[14, x.offset + c(7, 8, 9, 16)] <- 1 # baseline dissolution hazard predicted by degree distrib, cp, popgrowth
predictorMatrix[15, x.offset + c(7, 8, 9, 16)] <- 1 # age effect on dissolution hazard predicted by degree distrib, cp, popgrowth, HIV prev in older people (maybe?)
# NOTE: As it stands, each output statistic is predicted by ALL input and ALL other output statistics. That may not be a great idea, or even possible, if there is collinearity.
# Test dummy.MaC.incremental, and also dummy.MaC.incremental.parallel.mice
test.VEME2.MaC.incremental <- dummy.MaC.incremental.parallel.mice(targets.empirical = dummy.targets.empirical,
RMSD.tol.max = 0.95,
min.givetomice = 80, # 400
n.experiments = 320, # 1000
lls = c(1, 0.12, -0.3, 2.5, 0.1, 0.1, 20, 20, -0.8, 2, -0.35, -0.35, -3.6, -0.8, -0.16),
uls = c(1.2, 0.37, 0.3, 3.5, 0.4, 0.4, 66, 66, -0.25, 3.9, -0.1, -0.1, -1.4, -0.3, -0.001),
model = VEME.wrapper2.medians, # VEME.wrapper2,
strict.positive.params = c(4:8),
predictorMatrix = predictorMatrix,
maxit = 5,
maxwaves = 10,
n_cluster = 8) # 6
#(round(l1median(head(test.MaC.incremental$selected.experiments[[length(test.MaC.incremental$selected.experiments)]]), 1), 99)[5:8] - dummy.targets.empirical[1:4]) / dummy.targets.empirical[1:4]
#round(l1median(head(test.MaC.incremental$selected.experiments[[length(test.MaC.incremental$selected.experiments)]]), 1), 2)
#test.MaC.incremental$secondspassed
test.VEME2.MaC.incremental$secondspassed
test.VEME2.MaC.incremental$max.RMSD
test.VEME2.MaC.incremental$n.close.to.targets
head(test.VEME2.MaC.incremental$selected.experiments[[length(test.VEME2.MaC.incremental$selected.experiments)]])
save(dummy.targets.empirical, test.VEME2.MaC.incremental, file = "/Users/delvaw/Documents/MiceABC/test.VEME2.MaC.incremental.RData")
### Now we simulate for the 1 (or 5?) best fitting model(s)
inputs.calib <- as.numeric(test.VEME2.MaC.incremental$selected.experiments[[length(test.VEME2.MaC.incremental$selected.experiments)]][1, 1:15])
calib.dummy.master2 <- simpact.parallel(model = VEME.wrapper2,
actual.input.matrix = matrix(rep(inputs.calib, each = n.experiments), nrow = n.experiments),
seed_count = 0,
n_cluster = 8)
save(calib.dummy.master2, file = "/Users/delvaw/Documents/MiceABC/calib.dummy.master2.RData")
## AND THE OUTPUT LOOKS LIKE:
valid.inc.wide <- as.data.frame(valid.dummy.master2[, 17:31])
names(valid.inc.wide) <- 1:15 #paste0("incid.", 1:15)
valid.inc.wide$rep <- 1:nrow(valid.inc.wide)
valid.inc.wide$type <- "Master model"
calib.inc.wide <- as.data.frame(calib.dummy.master2[, 17:31])
names(calib.inc.wide) <- 1:15 #paste0("incid.", 1:15)
calib.inc.wide$rep <- (1 + nrow(valid.inc.wide)):(nrow(valid.inc.wide) + nrow(calib.inc.wide))
calib.inc.wide$type <- "Calibrated model"
inc.both <- rbind(valid.inc.wide, calib.inc.wide)
inc.both.long <- gather(inc.both, year, incidence, 1:15)
inc.both.long$year <- as.numeric(inc.both.long$year)
inc.both.long$smoother <- factor(paste0("smoother.", inc.both.long$type))
###
# Plotting the result
###
library(metafolio)
n.colours <- 2
cols <- gg_color_hue(n.colours, l=65)
darkcols <- gg_color_hue(n.colours, l=40)
incplot <- ggplot(filter(inc.both.long,
year >=10),
aes(year, incidence, group = rep, colour = type)) +
geom_line() +
facet_wrap(~ type) +
stat_smooth(se=FALSE,
method="loess",
span=1,
aes(year, incidence, group = type, colour=smoother)) +
xlab("Time (years)") +
ylab("HIV incidence (cases/person-year)") +
# scale_colour_hue(l = c(rep(65, 5), rep(10, 5))) +
scale_color_manual(values = c("Master model"=cols[1],
"Calibrated model"=cols[2],
"smoother.Master model"=darkcols[1],
"smoother.Calibrated model"=darkcols[2]),
guide = FALSE)
plot(incplot)
wdelva/MiceABC documentation built on May 4, 2019, 2:04 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#setwd("/Users/wimdelva/Documents/MiceABC/R")
#source("simpact.wrapper.R")
#source("VEME.wrapper.R")
#source("dummy.wrapper.R")
#source("simpact.parallel.R")
#source("MaC.incremental.R")
#source("dummy.MaC.incremental.R")
#source("/user/data/gent/vsc400/vsc40070/phylo/scripts/VEME.wrapper.R")
#source("/user/data/gent/vsc400/vsc40070/phylo/scripts/simpact.parallel.R")
#source("/user/data/gent/vsc400/vsc40070/phylo/scripts/dummy.MaC.incremental.R")
#source("/Users/delvaw/Documents/MiceABC/R/VEME.wrapper.R")
source("/Users/delvaw/Documents/MiceABC/R/VEME.wrapper2.R")
source("/Users/delvaw/Documents/MiceABC/R/VEME.wrapper2.medians.R")
source("/Users/delvaw/Documents/MiceABC/R/mice.wrapper.R")
source("/Users/delvaw/Documents/MiceABC/R/00-Functions.R")
source("/Users/delvaw/Documents/MiceABC/R/simpact.parallel.R")
source("/Users/delvaw/Documents/MiceABC/R/mice.parallel.R")
source("/Users/delvaw/Documents/MiceABC/R/dummy.MaC.incremental.R")
source("/Users/delvaw/Documents/MiceABC/R/dummy.MaC.incremental.parallel.mice.R")
library(dplyr)
library(MASS)
library(splines)
library(boot)
#library(haven)
library(ggplot2)
library(GGally)
library(fitdistrplus)
library(lmtest)
library(mclust)
#library(depth)
library(pcaPP)
#library(devtools)
#install_github("wdelva/RSimpactHelp")
library(RSimpactCyan)
library(RSimpactHelper)
library(lmtest)
library(mice)
#library(miceadds)
library(parallel)
library(randtoolbox)
library(EasyABC)
library(dplyr)
library(tidyr)
library(nlme)
library(lme4)
library(boot)
library(data.table)
library(parallel)
install.packages("adegenet", dep = TRUE) # Still to be done on vsc
library(adegenet) # Still to be done on vsc
library(ade4) # Still to be done on vsc
source("http://adegenet.r-forge.r-project.org/files/patches/auxil.R")
#install.packages("devtools")
#library(devtools)
#install_github("emvolz/treedater") # Still to be done on vsc
library(treedater) # Still to be done on vsc
dummy.input.vector <- c(1.1, 0.25, 0, 3, 0.23, 0.23, # what if 1.1 becomes 1.4
45, 45, #45, 45, # what if 45 becomes 60
-0.5, 2.8, -0.2, -0.2, -2.5, -0.52, -0.05)# c(1000, 2, 3, 4)
x.offset <- length(dummy.input.vector)
n.experiments <- 80
dummy.master2 <- simpact.parallel(model = VEME.wrapper2,
actual.input.matrix = matrix(rep(dummy.input.vector, each = n.experiments), nrow = n.experiments),
seed_count = 0,
n_cluster = 8)
save(valid.dummy.master2, file = "/Users/delvaw/Documents/MiceABC/valid.dummy.master2.RData")
#####
# The output of the master model
#####
#inc.master.vector only exists for the validation version of the master model
# inc.master.vector <- dummy.master2[, 17]
# #save(inc.master.vector, file = "/Users/delvaw/Documents/MiceABC/inc.master.vector.RData")
# #save(inc.master.vector.inflated, file = "/Users/delvaw/Documents/MiceABC/inc.master.vector.inflated.RData")
#
# hist(inc.master.vector[!is.na(inc.master.vector)], 14) # The distribution of HIV incidence after 10 years
# mean(inc.master.vector) # The mean HIV incidence after 10 years: 0.01527118
# median(inc.master.vector) # The median HIV incidence after 10 years: 0.01532902
# quantile(inc.master.vector, c(0.025, 0.975)) # 0.002462209 0.026417416
head(dummy.master2)
inc.master.vector <- dummy.master2[, 17]
new.infect.vector <- dummy.master2[, 18]
recent.ratio.vector <- dummy.master2[, 19]
mean(new.infect.vector, na.rm = TRUE) # The mean HIV incidence after 10 years
median(new.infect.vector, na.rm = TRUE) # The median HIV incidence after 10 years
quantile(new.infect.vector, c(0.025, 0.975), na.rm = TRUE)
# mean.br.len.vector <- dummy.master2[, 19]
hist(new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf], 20)
hist(recent.ratio.vector[!is.na(recent.ratio.vector) & recent.ratio.vector < Inf], 20)
plot(new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf],
inc.master.vector[!is.na(new.infect.vector) & new.infect.vector < Inf])
plot(recent.ratio.vector[inc.master.vector > 0.0142 & inc.master.vector < 0.0162],
inc.master.vector[inc.master.vector > 0.0142 & inc.master.vector < 0.0162])
cor.test(inc.master.vector[!is.na(new.infect.vector) & new.infect.vector < Inf],
new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf])
cor.test(inc.master.vector[!is.na(new.infect.vector) & new.infect.vector < Inf],
new.infect.vector[!is.na(new.infect.vector) & new.infect.vector < Inf])
#plot(mean.br.len.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf],
# inc.master.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf])
#cor.test(mean.br.len.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf],
# inc.master.vector[!is.na(mean.br.len.vector) & mean.br.len.vector < Inf])
#plot(mean.br.len.vector,
# new.infect.vector)
dummy.master2 <- dummy.master2 %>%
as.data.frame() %>%
dplyr::filter(complete.cases(.))
dummy.targets.empirical <- l1median(dummy.master2)
# dummy.targets.empirical (based on 400 repeats):
# 5.79782475 4.16132170 0.63770526 2.69246630 1.21359926 -1.56183331 0.25472579 2.09719359 0.02891330 0.13671790
# 0.06816378 0.34554390 0.31086783 0.37050293 0.41706055 1.01682475 7.07415631
dummy.targets.empirical <- c(5.79782475,
4.16132170,
0.63770526,
2.69246630,
1.21359926,
-1.56183331,
0.25472579,
2.09719359,
0.02891330,
0.13671790,
0.06816378,
0.34554390,
0.31086783,
0.37050293,
0.41706055,
1.01682475,
7.07415631)
# round(colMeans(dummy.master2), 3) # For interest sake, what are the marginal means?
predictorMatrix <- (1 - diag(1, length(c(dummy.input.vector, dummy.targets.empirical)))) # This is the default matrix.
# # Let's now modify the first 15 rows of this matrix, corresponding to the indicators of predictor variables for the input variables. In brackets the values for the master model.
predictorMatrix[1:x.offset, ] <- 0 # First we "empty" the relevant rows, then we refill them.
# We are currently not allowing input variables to be predicted by other predictor variables. Only via output variables. We could change this at a later stage.
predictorMatrix[1, x.offset + c(10, 11, 17)] <- 1 # relative susceptibility in young women is predicted by HIV prevalence in young men and women, and recent infections (~ incidence)
predictorMatrix[2, x.offset + 3] <- 1 # agescale predicted by slope
predictorMatrix[3, x.offset + c(1, 3, 6)] <- 1 # mean of the person-specific age gap preferences is predicted by slope, intercept and AAD
predictorMatrix[4, x.offset + c(2, 4, 5)] <- 1 # sd of the person-specific age gap preferences is predicted by SD, WSD, BSD
predictorMatrix[5, x.offset + c(7, 8, 9, 13, 16)] <- 1 # man gamma a predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.men, exp(growthrate)
predictorMatrix[6, x.offset + c(7, 8, 9, 12, 16)] <- 1 # woman gamma a predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.women, exp(growthrate)
predictorMatrix[7, x.offset + c(7, 8, 9, 13, 16, 17)] <- 1 # man gamma b predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.men, exp(growthrate), and recent infections (~ incidence)
predictorMatrix[8, x.offset + c(7, 8, 9, 12, 16, 17)] <- 1 # woman gamma b predicted by gamma shape.male, scale.male, pp.cp, hiv.prev.25.34.men, exp(growthrate), and recent infections (~ incidence)
predictorMatrix[9, x.offset + c(2, 4, 5, 7, 8, 14, 15, 16, 17)] <- 1 # formation.hazard.agegapry.gap_factor_x_exp is predicted by population growth, age gap variance, hiv prevalence, and recent infections (~ incidence)
predictorMatrix[10, x.offset + c(7, 8, 9, 12, 13, 16, 17)] <- 1 # baseline formation hazard predicted by HIV prevalence, cp, degree distrib. HIV prevalence, and recent infections (~ incidence)
predictorMatrix[11, x.offset + c(7, 8, 9, 12, 13, 16, 17)] <- 1 # numrel man penalty is predicted by degree distrib, cp, prev, popgrowth, and recent infections (~ incidence)
predictorMatrix[12, x.offset + c(7, 8, 9, 12, 13, 16, 17)] <- 1 # # numrel woman penalty is predicted by degree distrib, cp, prev, popgrowth, and recent infections (~ incidence)
predictorMatrix[13, x.offset + 16] <- 1 # conception.alpha_base is predicted by popgrowth
predictorMatrix[14, x.offset + c(7, 8, 9, 16)] <- 1 # baseline dissolution hazard predicted by degree distrib, cp, popgrowth
predictorMatrix[15, x.offset + c(7, 8, 9, 16)] <- 1 # age effect on dissolution hazard predicted by degree distrib, cp, popgrowth, HIV prev in older people (maybe?)
# NOTE: As it stands, each output statistic is predicted by ALL input and ALL other output statistics. That may not be a great idea, or even possible, if there is collinearity.
# Test dummy.MaC.incremental, and also dummy.MaC.incremental.parallel.mice
test.VEME2.MaC.incremental <- dummy.MaC.incremental.parallel.mice(targets.empirical = dummy.targets.empirical,
RMSD.tol.max = 0.95,
min.givetomice = 80, # 400
n.experiments = 320, # 1000
lls = c(1, 0.12, -0.3, 2.5, 0.1, 0.1, 20, 20, -0.8, 2, -0.35, -0.35, -3.6, -0.8, -0.16),
uls = c(1.2, 0.37, 0.3, 3.5, 0.4, 0.4, 66, 66, -0.25, 3.9, -0.1, -0.1, -1.4, -0.3, -0.001),
model = VEME.wrapper2.medians, # VEME.wrapper2,
strict.positive.params = c(4:8),
predictorMatrix = predictorMatrix,
maxit = 5,
maxwaves = 10,
n_cluster = 8) # 6
#(round(l1median(head(test.MaC.incremental$selected.experiments[[length(test.MaC.incremental$selected.experiments)]]), 1), 99)[5:8] - dummy.targets.empirical[1:4]) / dummy.targets.empirical[1:4]
#round(l1median(head(test.MaC.incremental$selected.experiments[[length(test.MaC.incremental$selected.experiments)]]), 1), 2)
#test.MaC.incremental$secondspassed
test.VEME2.MaC.incremental$secondspassed
test.VEME2.MaC.incremental$max.RMSD
test.VEME2.MaC.incremental$n.close.to.targets
head(test.VEME2.MaC.incremental$selected.experiments[[length(test.VEME2.MaC.incremental$selected.experiments)]])
save(dummy.targets.empirical, test.VEME2.MaC.incremental, file = "/Users/delvaw/Documents/MiceABC/test.VEME2.MaC.incremental.RData")
### Now we simulate for the 1 (or 5?) best fitting model(s)
inputs.calib <- as.numeric(test.VEME2.MaC.incremental$selected.experiments[[length(test.VEME2.MaC.incremental$selected.experiments)]][1, 1:15])
calib.dummy.master2 <- simpact.parallel(model = VEME.wrapper2,
actual.input.matrix = matrix(rep(inputs.calib, each = n.experiments), nrow = n.experiments),
seed_count = 0,
n_cluster = 8)
save(calib.dummy.master2, file = "/Users/delvaw/Documents/MiceABC/calib.dummy.master2.RData")
## AND THE OUTPUT LOOKS LIKE:
valid.inc.wide <- as.data.frame(valid.dummy.master2[, 17:31])
names(valid.inc.wide) <- 1:15 #paste0("incid.", 1:15)
valid.inc.wide$rep <- 1:nrow(valid.inc.wide)
valid.inc.wide$type <- "Master model"
calib.inc.wide <- as.data.frame(calib.dummy.master2[, 17:31])
names(calib.inc.wide) <- 1:15 #paste0("incid.", 1:15)
calib.inc.wide$rep <- (1 + nrow(valid.inc.wide)):(nrow(valid.inc.wide) + nrow(calib.inc.wide))
calib.inc.wide$type <- "Calibrated model"
inc.both <- rbind(valid.inc.wide, calib.inc.wide)
inc.both.long <- gather(inc.both, year, incidence, 1:15)
inc.both.long$year <- as.numeric(inc.both.long$year)
inc.both.long$smoother <- factor(paste0("smoother.", inc.both.long$type))
###
# Plotting the result
###
library(metafolio)
n.colours <- 2
cols <- gg_color_hue(n.colours, l=65)
darkcols <- gg_color_hue(n.colours, l=40)
incplot <- ggplot(filter(inc.both.long,
year >=10),
aes(year, incidence, group = rep, colour = type)) +
geom_line() +
facet_wrap(~ type) +
stat_smooth(se=FALSE,
method="loess",
span=1,
aes(year, incidence, group = type, colour=smoother)) +
xlab("Time (years)") +
ylab("HIV incidence (cases/person-year)") +
# scale_colour_hue(l = c(rep(65, 5), rep(10, 5))) +
scale_color_manual(values = c("Master model"=cols[1],
"Calibrated model"=cols[2],
"smoother.Master model"=darkcols[1],
"smoother.Calibrated model"=darkcols[2]),
guide = FALSE)
plot(incplot)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.