+ Näytä koodi- Piilota koodi
library(OpasnetUtils)
library(ggplot2)
objects.latest("Op_fi4305", code_name = "alusta") # [[Pneumokokkirokote]]
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]]
openv.setN(100)
## Read the annual IPD and carriage incidence data.
## The 0 entries in IPD and carriage data are replaced by small values.
serotypes<-c(
"19F", "23F", "6B", "14", "9V", "4", "18C", "1", "7",
"6A", "19A", "3", "8", "9N", "10", "11", "12", "15",
"16", "20", "22", "23A", "33", "35", "38", "6C", "Oth")
car_under5<-c(
156030, 156030, 126990, 41200, 22290, 12830, 10130, 10, 14180,
54940, 24320, 12160, 1350, 20940, 4050, 72270, 10, 33100,
3380, 1350, 12160, 3380, 680, 30400, 4050, 27470, 24320 )
car_over5<-c(
168100, 314800, 256700, 209800, 114100, 62500, 200700, 100, 100,
158800, 54900, 30800, 8800, 8800, 20800, 97700, 100, 100,
191900, 25200, 72500, 22000, 100, 71300, 100, 79400, 330100 )
ipd_under5<-c(
7.78, 7.88, 24.39, 20.76, 2.91, 2.91, 6.64, 0.31, 3.02,
3.94, 9.88, 1.25, 0.10, 0.83, 0.41, 0.42, 0.21, 1.98,
0.21, 0.01, 0.93, 0.10, 0.42, 0.31, 0.42, 0.01, 0.73 )
ipd_over5<-c(
28.51, 53.72, 29.53, 99.43, 43.07, 76.99, 24.39, 6.58, 46.88,
17.42, 20.54, 55.04, 11.21, 25.20, 6.28, 12.76, 13.89, 9.18,
4.73, 3.29, 29.03, 4.40, 5.64, 12.41, 1.43, 5.50, 11.20 )
## Combine the data into 2 matrices of dimension 27*2:
IPD<-cbind(ipd_under5, ipd_over5)
Car<-cbind(car_under5, car_over5)
## Row numbers corresponding to the 3 different PCV formulations
## in matrices IPD and Car. Note: there is no serotype 5 in our data.
pcv7rows<-seq(7); pcv10rows<-seq(9); pcv13rows<-seq(12)
## Example S1.2A: Calculate the predicted incidence of IPD for the non-vaccine
## types(NVTs) under PCV13. The predictions are calculated separately for the
## two age classes. These are the values reported on the bottom panel in
## Figure 2 (there given as per 100K incidences).
postvacc <-Vaccination(IPD,Car,VT_rows=pcv13rows,p=1,q=1)
cat("Incidence of invasive pneumococcal disease after vaccination program (p=1, q=1).\n")
oprint(cbind(serotypes,postvacc[[1]]))
cat("Number of carriers after vaccination program.\n")
oprint(cbind(serotypes,postvacc[[2]]))
## Example S1.2B: Decrease in IPD incidence after adding a single new serotype
## to PCV13 separately for the two age categories.
next_under5<-NextVT(IPD[,1],Car[,1], VT_rows=pcv13rows,p=1)
next_over5 <-NextVT(IPD[,2],Car[,2], VT_rows=pcv13rows,p=1)
cat("Next serotype to add to a vaccine (p=1).\n")
oprint(rbind(serotypes, next_under5, next_over5))
# Nämä taulukot kannattaisi transposata niin näyttäisivät siistimmiltä.
## Example S1.3A: The optimal sequence for under 5 year olds when replacement is 100%.
## The output shows the decreases in IPD incidence for each step,
## corresponding to Figure 5(C). The last serotype (row 27, the category "Other")
## is excluded from any vaccine composition but is taken into account as a
## replacing serotype at each stage.
opt<-OptimalSequence(IPD[,1],Car[,1],VT_rows=0,Excluded_rows=27,p=1.0,HowmanyAdded=20)
cat("Optimal sequence of serotypes to add to a vaccine (p=1, HowmanyAdded=20).\n")
oprint(rbind(serotypes[opt[1,]],opt[2,]))
## Example S1.3B: The optimal sequence for the whole population when
## replacement is 50% and the current composition includes the PCV7 serotypes.
opt<-OptimalSequence(IPD,Car, VT_rows=pcv7rows,Excluded_rows=length(serotypes),
p=0.5,HowmanyAdded=17)
cat("Optimal sequence of serotypes to add to a vaccine (p=0.5, HownamyAdded=17).\n")
oprint(rbind(serotypes[opt[1,]],opt[2,]))
###################################
## Read the annual IPD and carriage incidence data.
## The 0 entries in IPD and carriage data are replaced by small values.
IPD <- Ovariable("IPD", ddata = "Op_fi4305.pneumokokki_vaestossa")
IPD@data <- IPD@data[IPD@data$Observation == "Incidence" , colnames(IPD@data) != "Observation"]
Car <- Ovariable("Car", ddata = "Op_fi4305.pneumokokki_vaestossa")
Car@data <- Car@data[Car@data$Observation == "Carrier" , colnames(Car@data) != "Observation"]
servac <- Ovariable("servac", data = data.frame(
Vaccine = rep(c("Current", "New"), each = length(serotypes)),
Serotype = serotypes,
Result = as.numeric(c(
serotypes %in% c("19F", "23F", "6B", "14", "9V", "4", "18C", "1", "7","6A", "19A", "3"),
serotypes %in% servac_user
))
))
p <- Ovariable("p", data = data.frame(Result = p_user))
q <- EvalOutput(Ovariable("q", data = data.frame(Result = q_user)))
# EvalOutput must be used because q is mentioned twice in the code and there will otherwise be a merge mismatch.
# The true number of adult carriers may actually be larger than estimated. This adjusts for that.
Car <- Car * Ovariable("adjust", data = data.frame(Age = c("Under 5", "Over 5"), Result = c(1, adultcarriers)))
VacCar <- EvalOutput(VacCar)
VacIPD <- EvalOutput(VacIPD)
cat("Number of carriers\n")
oprint(summary(VacCar))
cat("Incidence of invasive pneumococcal disease.\n")
oprint(summary(VacIPD))
if("Iter" %in% colnames(VacCar@output)) N <- max(VacCar@output$Iter) else N <- 1
ggplot(VacCar@output, aes(x = Serotype, weight = result(VacCar) / N, fill = Vaccine)) + geom_bar(position = "dodge") + theme_gray(base_size = 24) +
labs(title = "Carriers", y = "Number of carriers in Finland")
ggplot(VacIPD@output, aes(x = Serotype, weight = result(VacIPD) / N, fill = Vaccine)) + geom_bar(position = "dodge") + theme_gray(base_size = 24) +
labs(title = "Incidence of invasive pneumococcal disease", y = "Number of cases / 100000 py")
ggplot(VacIPD@output, aes(x = Vaccine, weight = result(VacIPD) / N, fill = Age)) + geom_bar(position = "stack") + theme_gray(base_size = 24) +
labs(title = "Incidence of invasive pneumococcal disease", y = "Number of cases / 100000 py")
| |