Two-dimensional Monte Carlo
| Moderator:Jouni (see all) |
|
|
| Upload data
|
Question
How to perform two-dimensional Monte Carlo in Opasnet?
Answer
Rationale
See also
- This method is used by e.g. Health impact assessment
| Moderator:Jouni (see all) |
|
|
| Upload data
|
How to perform two-dimensional Monte Carlo in Opasnet?
#This is code Op_en7805/mc2d on page [[Two-dimensional Monte Carlo]]
library(OpasnetUtils)
# Funktio 2-ulotteisen Monte Carlon laskemiseen
# Esim dose on yksilökohtainen tieto mutta ERF on sama kaikille, joskin tuntematon.
# Siksi näitä kahta Iter-saraketta ei voi suoraan yhdistää. Tähän 2 ratkaisua:
## Ensin lasketaan henkilökohtaiset RR:t tms varsinaisella funktiolla,
### Sitten näistä arvotaan bootstrapillä uudet iteraatiot jotka aggregoidaan populaatiotasolle.
### Tässä versiossa on aina sama RR-Iter-kombinaatio.
## Ensin muutetaan yksilötason Iter Iter2:ksi, jolloin uusi indeksi paisuttaa ovariablen N-kertaiseksi.
### Sitten lasketaan ja lopuksi aggregoidaan populaatiotasolle.
### Tämä versio vaatii paljon enemmän muistia koska boostrap voidaan tehdä loopilla.
# Boostrap-versio:
# Paramter list. Note: this is not stored, you have to define it in the model code.
mc2dparam<- list(
N2 = 1000, # Number of iterations in the new Iter
run2d = TRUE, # Should the mc2d function be used or not?
newmarginals = c("Gender", "Ages", "Country"), # Names of columns that are non-marginals but should be sampled enough to become marginals
method = "bootstrap", # which method to use for 2D Monte Carlo? Currently bootsrap is the only option.
fun = mean # Function for aggregating the first Iter dimension.
)
mc2d <- function(ova) {
if(!exists("mc2dparam")) stop("Parameter list mc2dparam missing!\n")
require(reshape2)
marg <- setdiff(c(colnames(ova@output)[ova@marginal], mc2dparam$newmarginals), "Iter")
out <- aggregate(
result(ova),
by = ova@output[colnames(ova@output) %in% marg],
FUN = function(x) {
apply(
array(
as.numeric(sample(as.character(x), length(x)*mc2dparam$N2, replace=TRUE)),
#Numeric conversion is needed to prevent x from being interpreted as number of choices.
dim = c(length(x),mc2dparam$N)
),
MARGIN=2, FUN = mc2dparam$fun
)
}
)
temp <- melt(out[[length(out)]])
out[[length(out)]] <- 1:nrow(out)
colnames(temp) <- c("Nrow","Iter","Result")
out <- merge(out, temp, by.x = "x", by.y="Nrow")
out$x <- NULL
out <- Ovariable(
output = out,
marginal = colnames(out) %in% c(marg, "Iter")
)
return(out)
}
objects.store(mc2d)
cat("Function mc2d stored.\n")
|