Epidemiologinen malli

Laskenta

R-program code for the replcement model as in File S1 of Nurhonen M, Auranen K: Optimal serotype compositions for pneumococcal conjugate vaccination under serotype replacement [[1]

Esimerkkejä

VT: vaccine serotypes, i.e. pneumococcus serotypes that are found in a vaccine.
NVT: non-vaccine serotypes, i.e. pneumococcus serotypes that are not found in the vaccine.
Example model run

----#: . Yritin selvittää kuinka alla oleva koodi toimii. Huomasin että vaihtaessa arvoja ajoon, niin vaikuttaisi siltä että ajon tulokset eivät kuitenkaan muutu? --Jaakko (keskustelu) 5. kesäkuuta 2014 kello 14.13 (UTC) (type: truth; paradigms: science: comment) ----#: . Taulukot vastaavat näköjään aina lähtömallia, joka on PCv13, mutta pylväskuviossa muutettu malli on mukana. --Markku N. (keskustelu) 9. kesäkuuta 2014 kello 09.18 (UTC) (type: truth; paradigms: science: comment)

Funktioiden alustus

+ 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")

Data

Pneumokokin esiintyvyys suomalaisessa väestössä. Carrier: (oireettomien) kantajien lukumäärä, Incidence: invasiivisen pneumokokkitaudin ilmaantuvuus 100000 henkilövuotta kohti.

Pneumokokki väestössä(#, #/100000py)
Obs	Serotype	Age	Carrier	Incidence
1	19F	Under 5	156030	7.78
2	23F	Under 5	156030	7.88
3	6B	Under 5	126990	24.39
4	14	Under 5	41200	20.76
5	9V	Under 5	22290	2.91
6	4	Under 5	12830	2.91
7	18C	Under 5	10130	6.64
8	1	Under 5	10	0.31
9	7	Under 5	14180	3.02
10	6A	Under 5	54940	3.94
11	19A	Under 5	24320	9.88
12	3	Under 5	12160	1.25
13	8	Under 5	1350	0.1
14	9N	Under 5	20940	0.83
15	10	Under 5	4050	0.41
16	11	Under 5	72270	0.42
17	12	Under 5	10	0.21
18	15	Under 5	33100	1.98
19	16	Under 5	3380	0.21
20	20	Under 5	1350	0.01
21	22	Under 5	12160	0.93
22	23A	Under 5	3380	0.1
23	33	Under 5	680	0.42
24	35	Under 5	30400	0.31
25	38	Under 5	4050	0.42
26	6C	Under 5	27470	0.01
27	Oth	Under 5	24320	0.73
28	19F	Over 5	168100	28.51
29	23F	Over 5	314800	53.72
30	6B	Over 5	256700	29.53
31	14	Over 5	209800	99.43
32	9V	Over 5	114100	43.07
33	4	Over 5	62500	76.99
34	18C	Over 5	200700	24.39
35	1	Over 5	100	6.58
36	7	Over 5	100	46.88
37	6A	Over 5	158800	17.42
38	19A	Over 5	54900	20.54
39	3	Over 5	30800	55.04
40	8	Over 5	8800	11.21
41	9N	Over 5	8800	25.2
42	10	Over 5	20800	6.28
43	11	Over 5	97700	12.76
44	12	Over 5	100	13.89
45	15	Over 5	100	9.18
46	16	Over 5	191900	4.73
47	20	Over 5	25200	3.29
48	22	Over 5	72500	29.03
49	23A	Over 5	22000	4.4
50	33	Over 5	100	5.64
51	35	Over 5	71300	12.41
52	38	Over 5	100	1.43
53	6C	Over 5	79400	5.5
54	Oth	Over 5	330100	11.2

Epidemiologinen malli

Sisällys

Laskenta

Esimerkkejä

Funktioiden alustus

Data

Navigointivalikko