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]]
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 at 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
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 at
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 at 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