Assessment of the health impacts of H1N1 vaccination

This page is a assessment. The page identifier is Op_en4926
Moderator:Teemu R (see all)
This page is a stub. You may improve it into a full page, and then a rating bar will appear here.
Upload data Show results

Scope

What was the overall health impact of the H1N1 vaccination in Finland in 2009-2010? Given current knowledge, which was the better decision between vaccinating as happened versus vaccinating no-one?

Definition

Error creating thumbnail: Unable to save thumbnail to destination

Causal diagram.

Variables

• H1N1 vaccination coverage in Finland
• ERF of H1N1 vaccination on Narcolepsy
• A(H1N1)v immunity in the Finnish population
• Population of Finland
• Disability weights
• Life expectancy by age groups in Finland
• Probability of catching swine flu given subject is not immune
• Fraction of all cases represented by lab confirmed cases (which we have data on)
• Probability of death due to swine flu given a subject has swine flu and belongs to a risk group
• Fraction of population belonging to a risk group
• Length of swine flu

Indicators

• Narcolepsy in Finland
• AH1N1 cases in Finland

Formula

• Basic model
  • Uncertainties of ERF of vaccine on narcolepsy, fraction of all cases represented by lab confirmed cases and probability of catching swine flu are implemented.

+ Show code- Hide code # 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) # Function for aggregating indicators into DALYs # 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) outcome <- function(inarc = narc(), isf = sf(), inarcw = data.frame(Result=0.065),#op_baseGetData("opasnet_base", "Op_en2307")[,-c(1,2)], 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")[,-c(1)], ipop = op_baseGetData("opasnet_base", "Op_en2949")[,-c(1,2)], ivac_cov = op_baseGetData("opasnet_base", "Op_en4925")[,-c(1,2)], ibgi = data.frame(Result=1.4e-5), ...) { 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 iNERF 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 # 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; disabled sampling here since it would make different scenarios incomparable sf <- function(isfp = sfp(), ipop = op_baseGetData("opasnet_base", "Op_en2949")[,-c(1,2)], iimm = imm(), n = 1000, ilcf = tlcf, ...) { 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), hence variance is assumed to be variance^2 arbitrarily sfp <- function(dsf = op_baseGetData("opasnet_base", "Op_en4933")[,-c(1,2)], dpop = op_baseGetData("opasnet_base", "Op_en2949", series_id = 970)[,-c(1,2)], dimm = imm(), n = 1000, ...) { 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 a 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 sfpd <- function(dsf = op_baseGetData("opasnet_base", "Op_en4933")[,-c(1,2)], 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, ...) { 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) imm <- function(ibimm = op_baseGetData("opasnet_base", "Op_en4943", series_id = 994)[,-c(1,2)], ivac_cov = op_baseGetData("opasnet_base", "Op_en4925")[,-c(1,2)], ...) { 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) } # Calculating the outcome with different scenarios paralpha <- function(imean, ivar) imean * (imean * (1 - imean) / ivar - 1) parbeta <- function(imean, ivar) (1 - imean) * (imean * (1 - imean) / ivar - 1) library(OpasnetBaseUtils) library(MASS) library(ggplot2) # 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 # fast version below #vacscenario <- outcome(inarc = narc(iNERF = 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(1,rep(9.2, 3), rep(1, 13), 4)))) #novacscenario <- outcome(inarc = narc(iNERF = 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(1,rep(9.2, 3), rep(1, 13), 4)), # ivac_cov=data.frame(Result=0)), sf(iimm=imm(ivac_cov=data.frame(Result=0)))) #colnames(vacscenario)[ncol(vacscenario)] <- "vacscenario" #colnames(novacscenario)[ncol(novacscenario)] <- "novacscenario" #fout <- merge(vacscenario, novacscenario) #fout # slow ver below tvac_cov <- data.frame(op_baseGetData("opasnet_base", "Op_en4925")[,-c(1,2)], Scenario = "vacscenario") tvac_cov <- rbind(tvac_cov, data.frame(Age=tvac_cov[,"Age"], Result = 0, Scenario = "novacscenario")) test <- outcome(inarc = narc(ivac_cov = tvac_cov), sf(isfp = tsfp, iimm = imm(ivac_cov = tvac_cov))) 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) tapply(temp$Freq, temp[,c("Scenario")], mean) #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")[,-c(1)] 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 cor(temp2$Freq, temp2[,c("tNERF", "tlcf", "tsfp")], method = "spearman") ## EVPI mean(tapply(temp$Freq, temp$obs, min)) - min(tapply(temp$Freq, temp$Scenario, mean)) # not considering ## EVXPI for NERF, lcf and sfp 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)) - min(tapply(temp$Freq, temp[,c("Scenario")], mean)) 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)) - min(tapply(temp$Freq, temp[,c("Scenario")], mean)) 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)) - min(tapply(temp$Freq, temp[,c("Scenario")], mean)) #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) # What actually happened #defscenario <- outcome(inarc = op_baseGetData("opasnet_base", "Op_en4922")[,-1], isf = op_baseGetData("opasnet_base", "Op_en4933")[,-1]) # Narcolepsy erf = 1 #nerfscenario <- outcome(iNERF = data.frame(Result=1)) # No one was vaccinated #novacscenario <- outcome(ivac_cov = 0) # Other

Result

• From initial results it would appear like swine flu is more significant than narcolepsy in terms of DALYs.
  • Vaccinating as planned would result in approximately 2500 DALYs due to swine flu and narcolepsy combined.
  • Vaccinating no-one would result in approximately 4000 DALYs due to swine flu.