[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
Hello,
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) {
set.seed()
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)
return(dx)
}
i=1
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