[R] Convergence in Monte Carlo Simulation

[Michael Dewey]

I am not 100% clear what your code is doing as it gets a bit wangled as 
you posted in HTML but here are a couple of thoughts.

You need to set the seed outside any loops so it happens once and for all.

I would test after trycatch and keep a separate count of failures and 
successes as the failure to converge must be meaningful about the 
scientific question whatever that is. At the moment your count appears 
to be in the correct place to count successes.

Michael

On 14/06/2020 02:50, Phat Chau wrote:
> 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 using 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.
> 
> 

-- 
Michael
http://www.dewey.myzen.co.uk/home.html