a small improvement is to avoid computing the spectral decomposition
of `Sigma' (look at code of mvrnorm()) in each iteration, e.g.,
#set up a collector for subject means in A
a <- numeric(Ns)
#set up a collector for subject means in B
b <- numeric(Ns)
eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev <- eS$values
fact <- eS$vectors %*% diag(sqrt(pmax(ev, 0)), 2)
#start a monte carlo experiment
for (i in 1:mce) {
#generate correlated ideal means for each subject in each condition
X <- rnorm(2 * Ns)
dim(X) <- c(2, Ns)
sub.means <- t(fact %*% X)
#for each subject
for (s in 1:Ns) {
#generate some data for A and grab the mean
a[s] <- mean(rnorm(a.No[s], sub.means[s, 1], s.a.w[s]))
#generate some data for B and grab the mean
b[s] <- mean(rnorm(b.No[s], sub.means[s, 2], s.b.w[s]))
}
#store the observed correlation between A and B
sim.r[i] <- cor(a, b)
}
I hope it helps.
Best,
Dimitris
----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
Quoting Mike Lawrence <[EMAIL PROTECTED]>:
> Seems there were some line break issues when pasting the code, trying
> again with a different commenting style:
>
> #start a timer
> start = proc.time()[1]
>
> #set the true correllation
> rho = .5
>
> #set the number of Subjects
> Ns = 100
>
> #for each subject, set a number of observations in A
> a.No = 1:100
> #for each subject, set a number of observations in B
> b.No = 1:100
>
> #set the between Ss variance in condition A
> v.a = 1
> #set the between Ss variance in condition B
> v.b = 2
>
> #for each subject, set a standard deviation in A
> s.a.w = 1:100
> #for each subject, set a standard deviation in B
> s.b.w = 1:100
>
> #set the number of monte carlo experiments
> mce = 1e3
>
> #set up a collector for the simulated correlations
> sim.r = vector(length=mce)
>
> #define a covariance matrix for use in generating the correlated data
> Sigma=matrix(c(v.a,sqrt(v.a*v.b)*rho,sqrt(v.a*v.b)*rho,v.b),2,2)
>
> #set up a collector for subject means in A
> a = vector(length=Ns)
> #set up a collector for subject means in B
> b = vector(length=Ns)
>
> #start a monte carlo experiment
> for(i in 1:mce){
>
> #generate correlated ideal means for for each subject in each condition
> sub.means=mvrnorm(Ns,rep(0,2),Sigma)
>
> #for each subject
> for(s in 1:Ns){
>
> #generate some data for A and grab the mean
> a[s] = mean(rnorm(a.No[s],sub.means[s,1],s.a.w[s]))
> #generate some data for B and grab the mean
> b[s] = mean(rnorm(b.No[s],sub.means[s,2],s.b.w[s]))
>
> }
>
> #store the observed correlation between A and B
> sim.r[i] = cor(a,b)
>
> }
>
> #show the total time this took
> print(proc.time()[1]-start)
>
> ______________________________________________
> 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.
>
>
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
______________________________________________
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.
