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

From Testiwiki
Jump to: navigation, search
m (→‎Formula)
m (→‎Formula)
Line 117: Line 117:
 
# however if the probability has an uncertainty, only 1 sample from each distribution is sampled where
 
# 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
 
# samples of the probability will be used as inputs; disabled sampling here since it would make different scenarios incomparable
−
 
−
tlcf <- data.frame(Result = 0.1)#obs = 1:n, Result=rbeta(n, paralpha(0.1, 0.05^2), parbeta(0.1, 0.05^2)))
 
  
 
sf <- function(isfp = sfp(), ipop = op_baseGetData("opasnet_base", "Op_en2949")[,-c(1,2)], iimm = imm(), n = 1000,  
 
sf <- function(isfp = sfp(), ipop = op_baseGetData("opasnet_base", "Op_en2949")[,-c(1,2)], iimm = imm(), n = 1000,  
Line 149: Line 147:
 
# distribution. But we're interested in the uncertainty of the estimate... When n is suitably large a poisson distribution becomes  
 
# 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
 
# 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.
+
# 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",  
 
sfp <- function(dsf = op_baseGetData("opasnet_base", "Op_en4933")[,-c(1,2)], dpop = op_baseGetData("opasnet_base", "Op_en2949",  
Line 164: Line 165:
 
for (i in 1:nrow(out)) {
 
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]] /  
 
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),  
+
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]),  
+
parbeta(fitdistr(round(out$dsf[i]), "poisson")[[1]] / out$popnonimm[i], (fitdistr(round(out$dsf[i]),  
−
"poisson")[[2]]^2 / out$popnonimm[i]^2)))
+
"poisson")[[2]]^2 / out$popnonimm[i])^2)))
 
out2 <- rbind(out2, merge(out[i,], temp2, all.x = TRUE, all.y = FALSE))
 
out2 <- rbind(out2, merge(out[i,], temp2, all.x = TRUE, all.y = FALSE))
 
}
 
}
Line 191: Line 192:
 
out <- model.frame(I(dfrg * dsf / ilcf) ~., data = out)
 
out <- model.frame(I(dfrg * dsf / ilcf) ~., data = out)
 
colnames(out)[1] <- "Result"
 
colnames(out)[1] <- "Result"
−
out <- data.frame(Result=dmsf$dmsf/sum(out$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)
 
return(out)
 
}
 
}
Line 215: Line 218:
 
library(MASS)
 
library(MASS)
 
library(ggplot2)
 
library(ggplot2)
 +
 +
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)
 +
 +
#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
 
# Default, model based on current knowledge
Line 237: Line 248:
 
temp <- as.data.frame(as.table(tapply(test$Result, test[,c("obs","Scenario")], sum)))
 
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)
 
ggplot(temp, aes(x = Freq, fill = Scenario)) + geom_density(alpha = 0.2)
 +
tapply(temp$Freq, temp[,c("Scenario")], mean)
  
 
#vacscenario <- outcome()
 
#vacscenario <- outcome()

Revision as of 13:55, 5 April 2011



Scope

What was the overall health impact of the H1N1 vaccination, given that there is some evidence of it having caused narcolepsy.

Definition

Error creating thumbnail: Unable to save thumbnail to destination
Causal diagram.
Variables
Indicators

Formula

  • Basic model, works.
    • Uncertainties of ERF of vaccine on narcolepsy and probability of catching swine flu are implemented.

+ Show code

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.

Results

Conclusions

See also

References

  1. ↑ WHO life table estimates for 2008, Finland
Retrieved from "http://www.opasnet.org/testiwiki/index.php?title=Assessment_of_the_health_impacts_of_H1N1_vaccination&oldid=19971"