[R] Convergence in Monte Carlo Simulation

Phat Chau ph@t@ch@u @end|ng |rom m@||@utoronto@c@
Sun Jun 14 03:50:54 CEST 2020


I put together the following code and am curious about its correctness. My first question relates to the Monte Carlo simulations – the goal is to continue to iterate until I get n = 1000 simulations where the model successfully converges. I am wondering if I coded it correctly below with the while loop. Is the idea that the counter increments by one only if “model” does not return a string?

I would also like to know how I can create n = 1000 independent data sets. I think to do this, I would have to set a random number seed via set.seed() before the creation of each dataset. Where would I enter set.seed in the syntax below? Would it be in the function (as indicated in red)?

powercrosssw <- function(nclus, clsize) {


  cohortsw <- genData(nclus, id = "cluster")
  cohortsw <- addColumns(clusterDef, cohortsw)
  cohortswTm <- addPeriods(cohortsw, nPeriods = 8, idvars = "cluster", perName = "period")
  cohortstep <- trtStepWedge(cohortswTm, "cluster", nWaves = 4, lenWaves = 1, startPer = 1, grpName = "Ijt")

  pat <- genCluster(cohortswTm, cLevelVar = "timeID", numIndsVar = clsize, level1ID = "id")

  dx <- merge(pat[, .(cluster, period, id)], cohortstep, by = c("cluster", "period"))
  dx <- addColumns(patError, dx)

  setkey(dx, id, cluster, period)

  dx <- addColumns(outDef, dx)




while (i < 1000) {

  dx <- powercrosssw()

  #Fit multi-level model to simulated dataset
  model5 <- tryCatch(lme(y ~ factor(period) + factor(Ijt), data = dx, random = ~1|cluster, method = "REML"),
                     warning = function(w) { "warning" }

  if (! is.character(model5)) {

    coeff <- coef(summary(model5))["factor(Ijt)1", "Value"]
    pvalue <- coef(summary(model5))["factor(Ijt)1", "p-value"]
    error <- coef(summary(model5))["factor(Ijt)1", "Std.Error"]
    bresult <- c(bresult, coeff)
    presult <- c(presult, pvalue)
    eresult <- c(eresult, error)

    i <- i + 1

Thank you so much.

	[[alternative HTML version deleted]]

More information about the R-help mailing list