Difference between revisions of "Assessment of the health impacts of H1N1 vaccination"

# Model; original data inputs are disability weights (isfw, inarcw), population (ipop/dpop), (atm) ERF of vaccine on narcolepsy (NERF), 
# vaccination coverage (ivac_cov), base immunity in the population (ibimm) and observed number of sf cases (dsf)
# i in the beginning of the name of an object signifies 'input', d signifies 'data' and t signifies 'temporary'

# Function for aggregating indicators into DALYs. At the moment _narcolepsy cases_ and _swine flu cases_ are considered indicators, 
# however mortality due to swine flu should also be considered and indicator and calculated separately from this function. 
# Narcolepsy DALY weight taken as 0.065 (that of treated epilepsy)
# iEl is expectation of remaining life in an age group, value for "All" is arbitrary for now
# ilsf is the length of swine flu (e.g. how long was the subject absent frow work) in years
# ifrg is the fraction of population belonging to a risk group of dying from swine flu, default values are arbitrary
# iElrgc is the expected life correction for the risk group (people with heart conditions etc are likely to die earlier anyway)

outcome <- function(inarc = narc(), isf = sf(), inarcw = data.frame(Result=narcweight),#op_baseGetData("opasnet_base", "Op_en2307"), 
	isfw = data.frame(Result=0.5), iEl = data.frame(Age=c("0-4","5-9","10-14","15-19","20-24","25-29","30-34","35-39","40-44","45-49",
	"50-54","55-59","60-64","65-69","70-74","75-79","80+","All"), Result = c(79,75.1,70.1,65.1,60.3,55.5,50.7,45.9,41.1,36.5,31.9,27.6,
	23.5,19.5,15.7,12.1,9,50)), ilsf = data.frame(Result=5/365), isfpd = sfpd(), ifrg = data.frame(Age = c("0-4","5-9","10-14","15-19",
	"20-24","25-29","30-34","35-39","40-44","45-49","50-54","55-59","60-64","65-69","70-74","75-79","80+"), Result = c(0.1, rep(0.001,
	11), rep(0.5,5))), iElrgc = data.frame(Age = c("0-4","5-9","10-14","15-19","20-24","25-29","30-34","35-39","40-44","45-49","50-54",
	"55-59","60-64","65-69","70-74","75-79","80+"), Result = c(rep(1,12), rep(0.5,5))), ...) {
	
	names(inarc)[names(inarc)=="Result"] <- "inarc"
	names(inarcw)[names(inarcw)=="Result"] <- "inarcw"
	names(isf)[names(isf)=="Result"] <- "isf"
	names(isfw)[names(isfw)=="Result"] <- "isfw"
	names(iEl)[names(iEl)=="Result"] <- "iEl"
	names(ilsf)[names(ilsf)=="Result"] <- "ilsf"
	names(isfpd)[names(isfpd)=="Result"] <- "isfpd"
	names(ifrg)[names(ifrg)=="Result"] <- "ifrg"
	names(iElrgc)[names(iElrgc)=="Result"] <- "iElrgc"
	
	# Calculating narcolepsy DALYs
	onarc <- merge(inarc, inarcw)
	onarc <- merge(onarc, iEl)
	onarc <- model.frame(I(inarc * inarcw * iEl) ~., data = onarc)
	colnames(onarc)[1] <- "onarc"
	
	# Calculating swine flu DALYs from missing work
	osf <- merge(isf, isfw)
	osf <- merge(osf, ilsf)
	osf <- model.frame(I(isf * isfw * ilsf) ~., data = osf)
	colnames(osf)[1] <- "osf"
	
	#Calculating swine flu DALYs from deaths
	osfm <- merge(isf, isfpd)
	osfm <- merge(osfm, iEl)
	osfm <- merge(osfm, ifrg)
	osfm <- merge(osfm, iElrgc)
	osfm <- model.frame(I(isf * isfpd * iEl * ifrg * iElrgc) ~., data = osfm)
	colnames(osfm)[1] <- "osfm"
	
	# Aggregate (it is crucial that both variables are distinguished by the same age groups and other indices so that values
	# are not replicated to fill empty cells)
	out <- merge(onarc, osf)
	out <- merge(out, osfm)
	out <- reshape(out, idvar = names(out)[(names(out)%in%c("onarc","osf","osfm"))==FALSE], times = c("onarc","osf","osfm"), 
		timevar = "Outcome", varying = list(c("onarc","osf","osfm")), direction = "long")
	colnames(out)[ncol(out)] <- "Result"
	rownames(out) <- 1:nrow(out)
	
	return(out)
}

# Function for calculating narcolepsy cases, depends on erf of vac on narc, pop, vac_cov and background incidence (fraction)

narc <- function(iNERF = op_baseGetData("opasnet_base", "Op_en4923"), ipop = op_baseGetData("opasnet_base", "Op_en2949"), 
	ivac_cov = op_baseGetData("opasnet_base", "Op_en4925"), ibgi = data.frame(Result=1.4e-5), ...) {
	iNERF <- iNERF[, !colnames(iNERF)%in%c("id", "Result.Text")]
	ipop <- ipop[, !colnames(ipop)%in%c("id", "obs", "Result.Text")]
	ivac_cov <- ivac_cov[, !colnames(ivac_cov)%in%c("id", "obs", "Result.Text")]
	names(iNERF)[names(iNERF)=="Result"] <- "iNERF"
	names(ipop)[names(ipop)=="Result"] <- "ipop"	
	names(ivac_cov)[names(ivac_cov)=="Result"] <- "ivac_cov"
	names(ibgi)[names(ibgi)=="Result"] <- "ibgi"
	out <- merge(iNERF, ipop)
	out <- merge(out, ivac_cov)
	out <- merge(out, ibgi)
	out <- model.frame(I((iNERF - 1) * ipop * ivac_cov * ibgi) ~., data = out)
	colnames(out)[1] <- "Result"
	return(out)
}

# A function for calculating NERF might be added later

# Function for calculating swine flu cases
# ilcf is the fraction of actual swine flu cases represented by the lab confirmed cases
# (scrapped)assumed binomially distributed, if uncertainty of probability of getting swine flu is not 
# present, the point estimate will be used as such and n samples sampled from the same distribution;
# however if the probability has an uncertainty, only 1 sample from each distribution is sampled where
# samples of the probability will be used as inputs;(/scrapped) disabled sampling here since it would make different scenarios incomparable

sf <- function(isfp = sfp(), ipop = op_baseGetData("opasnet_base", "Op_en2949"), iimm = imm(), n = 1000, 
	ilcf = data.frame(Result = 0.2), ...) {
	ipop <- ipop[, !colnames(ipop)%in%c("id", "obs", "Result.Text")]
	names(isfp)[names(isfp)=="Result"] <- "isfp"
	names(ipop)[names(ipop)=="Result"] <- "ipop"	
	names(iimm)[names(iimm)=="Result"] <- "iimm"
	names(ilcf)[names(ilcf)=="Result"] <- "ilcf"
	out <- merge(isfp, ipop)
	out <- merge(out, iimm)
	out <- merge(out, ilcf)
	#Result <- NA
	#for (i in 1:nrow(out)) {
	#	if (sum(as.numeric(colnames(isfp)=="obs"))>0){
	#		Result[i] <- rbinom(1, round(out$ipop[i]*(1-out$iimm[i])), out$isfp[i]) / out$ilcf[i]
	#	} else {
	#		Result[((i-1)*n + 1):(i*n)] <- rbinom(n, round(out$ipop[i]*(1-out$iimm[i])), out$isfp[i]) / out$ilcf[i]
	#	}
	#}
	#out <- data.frame(out[,(colnames(out)%in%c("isfp","ipop","iimm","ilcf"))==FALSE], Result)
	out <- model.frame(I(isfp * ipop * (1 - iimm) / ilcf) ~., data = out)
	colnames(out)[1] <- "Result"
	return(out)
}

# Function for calculating probability of getting swine flu given subject is not vaccinated (simplified model)
# dsf, dpop and dimm are the actual data we have
# uncertainty: number of cases is binomially distributed, number of trials (n) is the number of people who were not immune to swine-flu
# outcome x is known, p is unknown. P can be easily calculated from x = p*n, when we assume that x is the expected value of the 
# distribution. But we're interested in the uncertainty of the estimate... When n is suitably large a poisson distribution becomes 
# an excellent approximation of the binomial distribution. So we can fit the poisson distribution using an R function (I didn't find
# one for binomial) from just one observation -> acquire value of x with standard error -> use as uncertainty of estimate; which produces
# rubbish because we only have a single observation, sd of the mean estimate of the poisson distribution appears to be (n*p)^0.5, so the 
# variance of the parameter is n * p, which divided by n gives p ... so the expected value is p and the variance is p, which causes problems
# with the beta distribution (shape parameters become negative), so I've used standard error^4 as the variance...
# Somebody with real knowledge about these things should check this sometime, I'm just playing around with some numbers to get more 
# uncertainty in the model.

sfp <- function(dsf = op_baseGetData("opasnet_base", "Op_en4933"), dpop = op_baseGetData("opasnet_base", "Op_en2949"), #series_id = 970),
	dimm = imm(), n = 1000, ...) {
	dsf <- dsf[, !colnames(dsf)%in%c("id", "obs", "Result.Text")]
	dpop <- dpop[, !colnames(dpop)%in%c("id", "obs", "Result.Text")]
	names(dsf)[names(dsf)=="Result"] <- "dsf"
	names(dpop)[names(dpop)=="Result"] <- "dpop"
	names(dimm)[names(dimm)=="Result"] <- "dimm"
	temp1 <- merge(dimm, dpop)
	temp1 <- model.frame(I(dpop * (1 - dimm)) ~., data = temp1)
	colnames(temp1)[1] <- "popnonimm"
	out <- merge(temp1, dsf)
	temp2 <- data.frame(obs=NA, Age=NA, fitp=NA)
	out2 <- data.frame()
	for (i in 1:nrow(out)) {
		temp2 <- data.frame(obs=1:n, Age=out$Age[i], fitp=rbeta(n, paralpha(fitdistr(round(out$dsf[i]), "poisson")[[1]] / 
			out$popnonimm[i], (fitdistr(round(out$dsf[i]), "poisson")[[2]]^2 / out$popnonimm[i])^2), 
			parbeta(fitdistr(round(out$dsf[i]), "poisson")[[1]] / out$popnonimm[i], (fitdistr(round(out$dsf[i]), 
			"poisson")[[2]]^2 / out$popnonimm[i])^2)))
		out2 <- rbind(out2, merge(out[i,], temp2, all.x = TRUE, all.y = FALSE))
	}
	out2 <- out2[,(colnames(out2)%in%c("dsf", "popnonimm"))==FALSE]
	colnames(out2)[ncol(out2)] <- "Result"
	return(out2)
}

# Function for calculating probability of mortality after catching swine flu given that the subject belongs to a risk group
# children of age < 1 and adults of 60+ afflicted by heart conditions etc are considered to belong to a risk group
# ifrg is the fraction belonging to the risk group in a given age group
# p(death|sf&rg)) assumed constant in all age groups
# dsf is the data we have on the number of swine flu cases, dmsf is the data on mortality caused by swine flu, dfrg is the data on 
# the fraction of population belonging to risk groups (arbitrary default values included)

sfpd <- function(dsf = op_baseGetData("opasnet_base", "Op_en4933"), dmsf = data.frame(Result=44), 
	dfrg = data.frame(Age = c("0-4 ","5-9","10-14","15-19","20-24","25-29","30-34","35-39","40-44","45-49",
	"50-54","55-59","60-64","65-69","70-74","75-79","80+"), Result = c(0.1, rep(0.001,11), rep(0.5,5))), ilcf = tlcf, ...) {
	dsf <- dsf[, !colnames(dsf)%in%c("id", "obs", "Result.Text")]
	names(dsf)[names(dsf)=="Result"] <- "dsf"
	names(dmsf)[names(dmsf)=="Result"] <- "dmsf"
	names(dfrg)[names(dfrg)=="Result"] <- "dfrg"
	names(ilcf)[names(ilcf)=="Result"] <- "ilcf"
	out <- merge(dsf, dfrg)
	out <- merge(out, ilcf)
	#out <- merge(out, dmsf)
	out <- model.frame(I(dfrg * dsf / ilcf) ~., data = out)
	colnames(out)[1] <- "Result"
	if (sum(as.numeric(colnames(out)=="obs"))>0) {#out <- as.data.frame(as.table(tapply(out$)))
		out <- data.frame(obs = levels(factor(out$obs)), Result = dmsf$dmsf/tapply(out$Result, out$obs, sum))
	} else {out <- data.frame(Result=dmsf$dmsf/sum(out$Freq))}
	return(out)
}

# Function for calculating total immunity to swine flu in the whole population
# P(imm) = P(vac or base_imm)
# ibimm is the input background immunity
# ivac_cov is the input vaccination coverage

imm <- function(ibimm = op_baseGetData("opasnet_base", "Op_en4943"), #series_id = 994), 
	ivac_cov = op_baseGetData("opasnet_base", "Op_en4925"), ...) {
	ibimm <- ibimm[, !colnames(ibimm)%in%c("id", "obs", "Result.Text")]
	ivac_cov <- ivac_cov[, !colnames(ivac_cov)%in%c("id", "obs", "Result.Text")]
	names(ibimm)[names(ibimm)=="Result"] <- "ibimm"
	names(ivac_cov)[names(ivac_cov)=="Result"] <- "ivac_cov"
	out <- merge(ibimm, ivac_cov)
	out <- model.frame(I(1 - (1 - ibimm) * (1 - ivac_cov)) ~., data = out)
	colnames(out)[1] <- "Result"
	return(out)
}

######################################################################3
## Calculating the outcome with different scenarios

# Libraries loaded and some basic functions defined

paralpha <- function(imean, ivar) imean * (imean * (1 - imean) / ivar - 1) # beta distribution parameters
parbeta <- function(imean, ivar) (1 - imean) * (imean * (1 - imean) / ivar - 1)
library(OpasnetUtils)
library(MASS)
library(ggplot2)
library(xtable)

# Common variable generation

#n <- 1000
tlcf <- data.frame(obs = 1:n, Result=rbeta(n, paralpha(0.2, 0.1^2), parbeta(0.2, 0.1^2)))
#ggplot(tlcf, aes(x=Result)) + geom_density()
#tlcf <- data.frame(Result = 0.1)

tsfp <- sfp()

#test2 <- data.frame(obs = 1:n, Result=rbeta(n, paralpha(0.2, 0.2^2+0.01), parbeta(0.2, 0.2^2+0.01)))
#ggplot(test2, aes(x=Result)) + geom_density()

# Default, model based on current knowledge

test <- op_baseGetData("opasnet_base", "Op_en4926") # Old run
test <- test[, !colnames(test)%in%c("id", "Result.Text")]
temp2 <- as.data.frame(as.table(apply(tapply(test$Result, test[,c("Age","Scenario","Outcome")], mean), c(2,3), sum)))
#ggplot(temp2, aes(x = Scenario, weight = Freq, fill = Outcome)) + geom_bar(position = "stack") # position = "dodge"
#head(temp2)
# For further analysis we need to use our recently generated lab confirmed fraction, this should be uploaded somewhere

temp <- op_baseGetData("opasnet_base", "Op_en4925") # New run. Op_en4925 = Vaccination coverage in Finland
temp <- temp[, !colnames(temp)%in%c("id", "obs", "Result.Text")]
tvac_cov <- data.frame(temp, Scenario = "Vaccinate all")
tvac_cov <- rbind(tvac_cov, data.frame(Age=tvac_cov[,"Age"], Result = 0, Scenario = "Vaccinate nobody"))
temp[2:4,2] <- 0
tvac_cov <- rbind(tvac_cov, data.frame(temp, Scenario = "Vaccinate all but 5-19 a"))

test <- outcome(inarc = narc(ivac_cov = tvac_cov), isf = sf(isfp = tsfp, iimm = imm(ivac_cov = tvac_cov), ilcf = tlcf))
test[test$Outcome=="onarc", "Outcome"] <- "Narcolepsy"
test[test$Outcome=="osf", "Outcome"] <- "Swine flu morbidity"
test[test$Outcome=="osfm", "Outcome"] <- "Swine flu mortality"
temp <- as.data.frame(as.table(tapply(test$Result, test[,c("obs","Scenario")], sum)))
ggplot(temp, aes(x = Freq, fill = Scenario)) + geom_density(alpha = 0.2)
temp4 <- as.data.frame(as.table(tapply(temp$Freq, temp[,c("Scenario")], mean)))
colnames(temp4) <- c("Decision option", "DALYs")
print(xtable(temp4), type = 'html')
ggplot(test, aes(x = Scenario, weight = Result, fill = Outcome)) + geom_bar(position = "stack") # position = "dodge"

#op_baseWrite("opasnet_base", test)

# Sensitivity and value of information analysis

#test <- op_baseGetData("opasnet_base", "Op_en4926")
temp2 <- as.data.frame(as.table(tapply(test$Result, test[,(colnames(test)%in%c("Result","Outcome"))==FALSE], sum)))
tNERF <- op_baseGetData("opasnet_base", "Op_en4923")
tNERF <- tNERF[, !colnames(tNERF)%in%c("id", "Result.Text")]
colnames(tNERF)[colnames(tNERF)=="Result"] <- "tNERF"
colnames(tlcf)[colnames(tlcf)=="Result"] <- "tlcf"#tlcf from above
colnames(tsfp)[colnames(tsfp)=="Result"] <- "tsfp"#tsfp from above
temp2 <- merge(temp2, tNERF)
temp2 <- merge(temp2, tlcf)
temp2 <- merge(temp2, tsfp)

temp2$NERFbin <- cut(temp2$tNERF, c(-Inf, 1, quantile(temp2$tNERF[temp2$tNERF>1], 1:9/10), Inf))
temp2$lcfbin <- cut(temp2$tlcf, c(-Inf, quantile(temp2$tlcf, 1:9/10), Inf))
temp2$sfpbin <- cut(temp2$tsfp, c(-Inf, quantile(temp2$tsfp, 1:9/10), Inf))

temp3 <- as.data.frame(as.table(tapply(temp2$Freq, temp2[,(colnames(temp2)%in%c("Freq","obs","tNERF","tlcf","tsfp",
	"lcfbin","sfpbin"))==FALSE], mean)))
colnames(temp3)[ncol(temp3)] <- "EogNERF"
temp2 <- merge(temp2, temp3)

temp3 <- as.data.frame(as.table(tapply(temp2$Freq, temp2[,(colnames(temp2)%in%c("Freq","obs","tNERF","tlcf","tsfp",
	"NERFbin","sfpbin","EogNERF"))==FALSE], mean)))
colnames(temp3)[ncol(temp3)] <- "Eoglcf"
temp2 <- merge(temp2, temp3)

temp3 <- as.data.frame(as.table(tapply(temp2$Freq, temp2[,(colnames(temp2)%in%c("Freq","obs","tNERF","tlcf","tsfp",
	"lcfbin","NERFbin","EogNERF","Eoglcf"))==FALSE], mean)))
colnames(temp3)[ncol(temp3)] <- "Eogsfp"
temp2 <- merge(temp2, temp3)

## Correlation

correlations <- as.data.frame(as.table(cor(temp2$Freq, temp2[,c("tNERF", "tlcf", "tsfp")], method = "spearman")))[, 2:3]
colnames(correlations) <- c("Variable", "Spearman correlation vs. outcome")
correlations$Variable <- c("Narcolepsy ERF", "Fraction of lab-confirmed cases", "P(swine flu|non-immune)")
print(xtable(correlations), type = 'html')

## EVPI

EVPI <- mean(tapply(temp$Freq, temp$obs, min)) - min(tapply(temp$Freq, temp$Scenario, mean)) # not considering by age group

EVPI.agegroup <- mean(apply(tapply(temp2$Freq, temp2[,c("obs","Age")], min), 1, sum)) - sum(apply(tapply(temp2$Freq, temp2[,c("Scenario","Age")], mean), 2, min)) # Not used in the table because I am not quite sure what this means

## EVXPI for NERF, lcf and sfp

EVXPI.NERF <- mean(apply(tapply(temp2$EogNERF, temp2[,(colnames(temp2)%in%c("Freq","Scenario","tNERF","tlcf","tsfp","NERFbin","lcfbin","sfpbin",
	"EogNERF","Eoglcf","Eogsfp"))==FALSE], min), 2, sum)) - sum(apply(tapply(temp2$Freq, temp2[,c("Scenario","Age")], mean), 2, min))

EVXPI.lcf <- mean(apply(tapply(temp2$Eoglcf, temp2[,(colnames(temp2)%in%c("Freq","Scenario","tNERF","tlcf","tsfp","NERFbin","lcfbin","sfpbin",
	"EogNERF","Eoglcf","Eogsfp"))==FALSE], min), 2, sum)) - sum(apply(tapply(temp2$Freq, temp2[,c("Scenario","Age")], mean), 2, min))

EVXPI.sfp <- mean(apply(tapply(temp2$Eogsfp, temp2[,(colnames(temp2)%in%c("Freq","Scenario","tNERF","tlcf","tsfp","NERFbin","lcfbin","sfpbin",
	"EogNERF","Eoglcf","Eogsfp"))==FALSE], min), 2, sum)) - sum(apply(tapply(temp2$Freq, temp2[,c("Scenario","Age")], mean), 2, min))

VOI <- data.frame(VOI = c("Total VOI (EVPI)", "Narcolepsy ERF", "Fraction of lab-confirmed cases", "P(swine flu|non-immune)"), value = c(EVPI, EVXPI.NERF, EVXPI.lcf, EVXPI.sfp))
colnames(VOI)[2] <- "Value in DALYs"
print(xtable(VOI), type = 'html')

#vacscenario <- outcome()
#novacscenario <- outcome(inarc = narc(ivac_cov=data.frame(Result=0)), sf(iimm=imm(ivac_cov=data.frame(Result=0))))
#temp <- as.data.frame(as.table(tapply(vacscenario$Result, vacscenario[,c("obs","Outcome")], sum)))
#tapply(vacscenario$Result, vacscenario[,c("Age","Outcome")], mean)
#tapply(temp$Freq, temp$Outcome, mean)
#ggplot(temp, aes(x=Freq, fill=Outcome)) + geom_density(alpha = 0.2)
#colnames(vacscenario)[ncol(vacscenario)] <- "vacscenario"
#colnames(novacscenario)[ncol(novacscenario)] <- "novacscenario"
#fout <- merge(vacscenario, novacscenario)
#fout <- reshape(fout, idvar = names(fout)[(names(fout)%in%c("vacscenario","novacscenario"))==FALSE], times = c("vacscenario","novacscenario"), 
#	timevar = "Scenario", varying = list(c("vacscenario","novacscenario")), direction = "long")
#colnames(fout)[ncol(fout)] <- "Result"
#temp <- as.data.frame(as.table(tapply(fout$Result, fout[,c("obs","Scenario")], sum)))
#ggplot(temp, aes(x = Freq, fill = Scenario)) + geom_density(alpha = 0.2)