Hello,
I have a simple (?) simulation problem.
I'm doing a simulation with logistic model and I want to reapet it 600
times.
The simulation looks like this:
z <- 0
x <- 0
y <- 0
aps <- 0
tiss <- 0
for (i in 1:500){
z[i] <- rbinom(1, 1, .6)
x[i] <- rbinom(1, 1, .95)
y[i] <- z[i]*x[i]
if (y[i]==1) aps[i] <- rnorm(1,mean=13.4, sd=7.09) else aps[i] <-
rnorm(1,mean=12.67, sd=6.82)
if (y[i]==1) tiss[i] <- rnorm(1,mean=20.731,sd=9.751) else tiss[i] <-
rnorm(1,mean=18.531,sd=9.499)
}
v <- data.frame(y, aps, tiss)
log_v <- glm(y~., family=binomial, data=v)
summary(log_v)
I want to do a repetition of this 600 times (I want to have 600 logistic
models), and see all the coefficients of the covariates aps & tiss.
Thanks in advance,
Sigalit.
[[alternative HTML version deleted]]
Dear Sigalit, On Nov 12, 2007 2:18 PM, sigalit mangut-leiba <smangut at gmail.com> wrote:> Hello, > I have a simple (?) simulation problem. > I'm doing a simulation with logistic model and I want to reapet it 600 > times. > The simulation looks like this: > > z <- 0 > x <- 0 > y <- 0 > aps <- 0 > tiss <- 0 > for (i in 1:500){ > z[i] <- rbinom(1, 1, .6) > x[i] <- rbinom(1, 1, .95) > y[i] <- z[i]*x[i] > if (y[i]==1) aps[i] <- rnorm(1,mean=13.4, sd=7.09) else aps[i] <- > rnorm(1,mean=12.67, sd=6.82) > if (y[i]==1) tiss[i] <- rnorm(1,mean=20.731,sd=9.751) else tiss[i] <- > rnorm(1,mean=18.531,sd=9.499) > } > v <- data.frame(y, aps, tiss) > log_v <- glm(y~., family=binomial, data=v) > summary(log_v) > > I want to do a repetition of this 600 times (I want to have 600 logistic > models), and see all the coefficients of the covariates aps & tiss. > Thanks in advance, > Sigalit. >I won't comment on your general approach to simulation but does this answer your question? ### Create an empty matrix to hold the results as they ### are created Results <- matrix(NA, ncol = 3, nrow = 600) for(i in 1:600) { generate your data Results[i, ] <- coef(glm(y~., family=binomial, data=v)) } summary(Results) Stephen -- Rochester, Minn. USA
Thank you, Sigalit. On 11/13/07, Moshe Olshansky <m_olshansky@yahoo.com> wrote:> > Hi, > > Look at names(log_v) in your notation to see what you > really need. > Then you could do something like: > > results <- list(600) > for (ij in 1:600) { > do what you did > results[[ij]] <- log_v > } > > So now you have a list of 600 results and can do > anything you need (look at lapply/sapply). > > Regards, > > Moshe. > > --- sigalit mangut-leiba <smangut@gmail.com> wrote: > > > Hello, > > I have a simple (?) simulation problem. > > I'm doing a simulation with logistic model and I > > want to reapet it 600 > > times. > > The simulation looks like this: > > > > z <- 0 > > x <- 0 > > y <- 0 > > aps <- 0 > > tiss <- 0 > > for (i in 1:500){ > > z[i] <- rbinom(1, 1, .6) > > x[i] <- rbinom(1, 1, .95) > > y[i] <- z[i]*x[i] > > if (y[i]==1) aps[i] <- rnorm(1,mean=13.4, sd=7.09) > > else aps[i] <- > > rnorm(1,mean=12.67, sd=6.82) > > if (y[i]==1) tiss[i] <- > > rnorm(1,mean=20.731,sd=9.751) else tiss[i] <- > > rnorm(1,mean=18.531,sd=9.499) > > } > > v <- data.frame(y, aps, tiss) > > log_v <- glm(y~., family=binomial, data=v) > > summary(log_v) > > > > I want to do a repetition of this 600 times (I want > > to have 600 logistic > > models), and see all the coefficients of the > > covariates aps & tiss. > > Thanks in advance, > > Sigalit. > > > > [[alternative HTML version deleted]] > > > > ______________________________________________ > > R-help@r-project.org mailing list > > 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. > > > >[[alternative HTML version deleted]]
And what is your question? Julian sigalit mangut-leiba wrote:> Hello, > I have a simple (?) simulation problem. > I'm doing a simulation with logistic model and I want to reapet it 600 > times. > The simulation looks like this: > > z <- 0 > x <- 0 > y <- 0 > aps <- 0 > tiss <- 0 > for (i in 1:500){ > z[i] <- rbinom(1, 1, .6) > x[i] <- rbinom(1, 1, .95) > y[i] <- z[i]*x[i] > if (y[i]==1) aps[i] <- rnorm(1,mean=13.4, sd=7.09) else aps[i] <- > rnorm(1,mean=12.67, sd=6.82) > if (y[i]==1) tiss[i] <- rnorm(1,mean=20.731,sd=9.751) else tiss[i] <- > rnorm(1,mean=18.531,sd=9.499) > } > v <- data.frame(y, aps, tiss) > log_v <- glm(y~., family=binomial, data=v) > summary(log_v) > > I want to do a repetition of this 600 times (I want to have 600 logistic > models), and see all the coefficients of the covariates aps & tiss. > Thanks in advance, > Sigalit. > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at r-project.org mailing list > 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.
Thank you, I changed that and it's much more efficient. Sigalit. On 11/13/07, Phil Spector <spector@stat.berkeley.edu> wrote:> > You can use the replicate function to do your simulation. First, > put the code to do one repetition in a function: > > dosim0 = function(n=500){ > x <- 0 > y <- 0 > z <- 0 > aps <- 0 > tiss <- 0 > for (i in 1:n){ > z[i] <- rbinom(1, 1, .6) > x[i] <- rbinom(1, 1, .95) > y[i] <- z[i]*x[i] > if (y[i]==1) aps[i] <- rnorm(1,mean=13.4, sd=7.09) else aps[i] <- > rnorm(1,mean=12.67, sd=6.82) > if (y[i]==1) tiss[i] <- rnorm(1,mean=20.731,sd=9.751) else tiss[i] > <- rnorm(1,mean=18.531,sd=9.499) > } > v <- data.frame(y, aps, tiss) > glm(y~., family=binomial, data=v)$coef > } > > > Now you can run > > results = t(replicate(1000,dosim0())) > > to get a 1000x3 matrix with the coefficients from each simulation. > > Before you actually do the simulations, you might want to rewrite > your function so that it's not so inefficient. Functions in R are > vectorized, and you should take advantage of that whenever possible. > Here's a vectorized version of your simulation function: > > dosim = function(n=500){ > z = rbinom(n,1,.6) > x = rbinom(n,1,.95) > y = z * x > aps = ifelse(y == 1,rnorm(n,mean=13.4,sd=7.09),rnorm(n,mean=12.67, sd> 6.82)) > tiss = ifelse(y == 1,rnorm(n,mean=20.731,sd=9.751),rnorm(n,mean=18.531 > ,sd=9.499)) > glm(y~.,family=binomial,data=data.frame(y,aps,tiss))$coef > } > > Now, you can run > > results = t(replicate(1000,dosim())) > > to get a similar result. Notice the difference in time for the > two functions: > > > system.time(results <- t(replicate(1000,dosim0()))) > user system elapsed > 55.123 0.004 55.239 > > system.time(results <- t(replicate(1000,dosim()))) > user system elapsed > 20.297 0.000 20.295 > > - Phil Spector > Statistical Computing Facility > Department of Statistics > UC Berkeley > spector@stat.berkeley.edu > >[[alternative HTML version deleted]]