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