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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
