Difference between revisions of "Two-dimensional Monte Carlo"

#This is code Op_en####/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:

mc2dparam<- list(
  N2 = 1000,
  run2d = TRUE,
  newmarginals = c("Gender", "Ages", "Country"),
  method = "bootstrap",
  fun = mean
)

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.latest(mc2d)
cat("Function mc2d stored.\n")