[R] Convergence in Monte Carlo Simulation

[Michael Dewey]

Dear Edward

Every time you call your function powercrosssw() it resets the seed so 
you must be calling it multiple times in some way.

Michael

On 14/06/2020 13:57, Phat Chau wrote:
> Thank you Michael.
> 
> I will clarify some more. The function in the first part of the code that I posted generates the simulated dataset for a cluster randomized trial from the simstudy package.
> 
> I am not quite clear what you mean by placing it outside the loop. So the goal here is to create n = 1000 independent datasets with different (randomly drawn values from the specified normal distributions not shown) for all of the parameters. What I have tried to do is place the seed at the very top of all my code in the past, but what that does is it leads to the creation of a single dataset that gets repeated over and over n = 1000 times. Hence, there ends up being no variability in the data (and power estimates from the p-values given the stated and required power).
> 
> Regarding the counter, is it correct in this instance that the loop will continue until n = 1000 iterations have successfully converged? I am not so concerned with counting failures.
> 
> Thank you.
> Edward
> 
> On 2020-06-14, 6:46 AM, "Michael Dewey" <lists using dewey.myzen.co.uk> wrote:
> 
>      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
>      
> 
> 
> 

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