+ Näytä koodi- Piilota koodi
library(OpasnetUtils)
#S1.4. The R-functions
###############################################################################
##
## R code for the core methods introduced in
## Markku Nurhonen and Kari Auranen:
## "Optimal serotype compositions for pneumococcal conjugate
## vaccination under serotype replacement",
## PLoS Computational Biology, 2014.
##
###############################################################################
## List of arguments common to most functions:
##
## IPD = matrix of IPD incidences by age class (columns) and serotype (rows)
## Car = corresponding matrix of carriage incidences
## VT_rows = vector of the row numbers in matrices IPD and Car
## corresponding to vaccine types (VT_rows=0 for no vaccination)
## p = proportion of lost VT carriage which is replaced by NVT carriage
## q = proportion of VT carriage lost either due to elimination or replacement
##
## This code includes 4 functions:
## Vaccination, NextVT, OptimalSequence and OptimalVacc.
##
Vaccination<-function(IPD,Car,VT_rows,p,q) {
##
## Result:
## A list of 2 matrices: IPD and carriage incidences
## after vaccination (corresponding to matrices IPD and Car).
## [Markku Nurhonen 2013]
##
if (VT_rows[1]>0) {
IPD<-as.matrix(IPD); Car<-as.matrix(Car)
# Post vaccination carriage incidences
Car_Total<-t(matrix(apply(Car,2,sum),dim(Car)[2],dim(Car)[1]))
Car2<-Car; Car2[VT_rows,]<-0
Car_NVT<-t(matrix(apply(Car2,2,sum),dim(Car2)[2],dim(Car2)[1]))
Car_VT<-Car_Total-Car_NVT
CarNew<-q*(1+p*Car_VT/Car_NVT)*Car2+(1-q)*Car
# Post vaccination IPD incidences
NVT_rows<-seq(dim(IPD)[1])[-1*VT_rows]
# CCR=Case-to-carrier ratios
CCR<-IPD/Car ; IPDNew<-0*IPD
# Apply the equation appearing above
# equation (1) in text for each serotype.
# First term applies to NVTs.
IPDNew[VT_rows,]<-(1-q)*IPD[VT_rows,]
# Second term applies to NVTs.
IPDNew[NVT_rows,]<-((Car_NVT+p*q*Car_VT)*(Car/Car_NVT)*CCR)[NVT_rows,]
}
else {
IPDNew<-IPD; CarNew<-Car
}
list(IPDNew,CarNew)
}
NextVT<-function(IPD,Car,VT_rows,p) {
##
## Result:
## A vector of decreases in IPD due to adding a serotype
## to the vaccine. If VT_rows=0, initially no vaccination.
## For row indexes incuded in VT_rows, the result is 0.
## [Markku Nurhonen 2013]
##
IPD<-as.matrix(IPD); Car<-as.matrix(Car)
## VaccMat = IPD and Car matrices after vaccination
VaccMat<-Vaccination(IPD,Car,VT_rows,p,1)
IPD<-VaccMat[[1]]; Car<-VaccMat[[2]]
## Total_IPD,Total_Car = Matrices corresponding to
## overall IPD and carriage in each age class.
Total_IPD<-t(matrix(apply(IPD,2,sum),dim(IPD)[2],dim(IPD)[1]))
Total_Car<-t(matrix(apply(Car,2,sum),dim(Car)[2],dim(Car)[1]))
## Effect = decrease in IPD when one serotype is added to the vaccine.
## See equation (3) in text.
Effect<-(Total_IPD-IPD)*((IPD/(Total_IPD-IPD))-(p*Car/(Total_Car-Car)))
## Special case when only one NVT remains.
IPD_nonzero<-which(apply(IPD,1,sum)!=0)
if (length(IPD_nonzero)==1) {Effect[IPD_nonzero,]<-IPD[IPD_nonzero,]}
## Result is obtained after summation over age classes.
apply(Effect,1,sum)
}
OptimalSequence<-function(IPD,Car,VT_rows,Excluded_rows,p,HowmanyAdded) {
##
## Starting from VTs indicated by the vector VT_rows
## (VT_rows=0, for no vaccination) sequentially add new VTs
## to the vaccine composition s.t. at each step the optimal
## serotype (corresponding to largest decrease in IPD) is added.
##
## Excluded_rows = Vector of indexes of the rows in matrices
## IPD and Car corresponding to serotypes that are not to
## be included in a vaccine composition, e.g. a row
## corresponding to a group of serotypes labelled "Other".
## Enter Excluded_rows=0 for no excluded serotypes.
## HowmanyAdded = number of VTs to be added.
##
## Result:
## Matrix of dimension 2*HowmanyAdded with 1st row indicating
## the row numbers of added serotypes in the order they appear
## in the sequence. The 2nd row lists the decreases in IPD
## due to addition of each type. [Markku Nurhonen 2013]
##
IPD<-as.matrix(IPD); Car<-as.matrix(Car)
## First check the maximum possible number of added VTs.
VT_howmany<-length(VT_rows)
if (VT_rows[1]==0) {VT_howmany<-0}
Excluded_howmany<-length(Excluded_rows)
if (Excluded_rows[1]==0) {Excluded_howmany<-0}
HowmanyAdded<-min(HowmanyAdded,dim(IPD)[1]-(VT_howmany+Excluded_howmany))
BestVTs<-BestEffects<-rep(0,HowmanyAdded)
## Sequential procedure: at each step find the best additional VT.
for (i in 1:HowmanyAdded) {
## Effects = Decrease in IPD after addition of each serotype
Effects<-NextVT(IPD,Car,VT_rows,p)
## Set Effects for VTs and excluded types equal to small values
## so that none of these will be selected as the next VT.
minvalue<- -2*max(abs(Effects))
if (Excluded_howmany>0) {Effects[Excluded_rows]<-minvalue}
if (VT_rows[1]>0) {Effects[VT_rows]<-minvalue}
## BestVTs[i] = Index of serotype with maximum decrease in IPD.
BestVTs[i]<-order(-1*Effects)[1]
## BestEffects[i] = Decrese in IPD due to addition of BestVTs[i]
## to the vaccine.
BestEffects[i]<-Effects[BestVTs[i]]
VT_rows<-c(VT_rows,BestVTs[i])
if (VT_rows[1]==0) {VT_rows<-VT_rows[-1]}
VaccMat<-Vaccination(IPD,Car,VT_rows,p,1)
IPD<-VaccMat[[1]]; Car<-VaccMat[[2]]
}
t(matrix(c(BestVTs,BestEffects),HowmanyAdded,2))
}
OptimalVacc<-function(IPD,Car,VT_rows,p,q,HowmanyAdded) {
##
## Result:
## A list of 3 elements: (1) Row numbers of serotypes in the optimal
## vaccine composition (2)-(3) IPD and carriage incidences
## by serotype and age class corresponding to the optimal
## vaccine formed using the sequential procedure in the
## function OptimalSequence. [Markku Nurhonen 2013]
##
Additional_VTs<-OptimalSequence(IPD,Car,VT_rows,p,HowmanyAdded)[1,]
All_VTs<-c(VT_rows,Additional_VTs)
if (All_VTs[1]==0) All_VTs<-All_VTs[-1]
VaccMat<-Vaccination(IPD,Car,All_VTs,p,q)
list(All_VTs,VaccMat[[1]],VaccMat[[2]])
}
VacCar <- Ovariable("VacCar",
dependencies = data.frame(Name = c(
"IPD", # incidence of pneumococcus disease
"Car", # number of carriers of pneumococcus
"servac", # ovariable of serotypes in vaccine (1 for serotypes in a vaccine, otherwise result is 0)
"p", # proportion of eliminated VT carriage that is replaced by NVT carriage
"q" # proportion of of VT carriage eliminated by vaccine
)),
formula = function(...) {
## Result:
## An ovariable of carriage incidences
## after vaccination (corresponding to Car).
## [Markku Nurhonen 2013, Jouni Tuomisto 2014]
# Post vaccination carriage incidences
# Sum over serotypes and drop extra columns
Car_Total<- unkeep(oapply(Car, cols = "Serotype", FUN = sum) * 1, prevresults = TRUE)
# Car2 is a temporary ovariable with NVT carriers only
Car2 <- unkeep(Car * (1 - servac), prevresults = TRUE) # Take only NVT carriers
Car_NVT <- oapply(Car2, cols = "Serotype", FUN = sum) # Carriers of serotypes not in vaccine (NVT)
Car_VT <- Car_Total - Car_NVT # Carriers of vaccine serotypes
CarNew <- q * (1 + p * Car_VT / Car_NVT) * Car2 + (1 - q) * Car
return(CarNew)
}
)
VacIPD <- Ovariable("VacIPD",
dependencies = data.frame(Name = c(
"IPD", # incidence of pneumococcus disease
"Car", # number of carriers of pneumococcus
"servac", # ovariable of serotypes in vaccine (1 for serotypes in a vaccine, otherwise result is 0)
"p", # proportion of eliminated VT carriage that is replaced by NVT carriage
"q" # proportion of of VT carriage eliminated by vaccine
)),
formula = function(...) {
## Result:
## An ovariable of IPD incidence
## after vaccination (corresponding to ovariable IPD).
## [Markku Nurhonen 2013, Jouni Tuomisto 2014]
# Post vaccination carriage incidences (same code as in VacCar)
Car_Total <- unkeep(oapply(Car, cols = "Serotype", FUN = sum) * 1, prevresults = TRUE) # Sums over serotypes
Car2 <- unkeep(Car * (1 - servac), prevresults = TRUE)
Car_NVT <- oapply(Car2, cols = "Serotype", FUN = sum) # Carriers of serotypes not in vaccine (NVT)
Car_VT <- Car_Total - Car_NVT # Carriers of vaccine serotypes
CarNew <- q * (1 + p * Car_VT / Car_NVT) * Car2 + (1 - q) * Car
# Post vaccination IPD incidences
# CCR=Case-to-carrier ratios
CCR <- IPD / Car
# Apply the equation appearing above
# equation (1) in text for each serotype.
# First term applies to VTs.
IPDNewVT <- (1 - q) * IPD * servac
# Second term applies to NVTs.
IPDNewNVT <- (Car_NVT + p * q * Car_VT) * (Car / Car_NVT) * CCR * (1 - servac)
IPDNew <- IPDNewVT + IPDNewNVT
return(IPDNew)
}
)
objects.store(Vaccination, NextVT, OptimalSequence, OptimalVacc, VacCar, VacIPD)
cat("Funktiot Vaccination, NextVT, OptimalSequence, OptimalVacc sekä ovariablet VacCar, VacIPD tallennettu. \n")
| |
