Dear List:

First, many thanks to those who offered assistance while I constructed
code for the simulation. I think I now have code that resolves most of
the issues I encountered with memory.

While the code works perfectly for smallish datasets with small sample
sizes, it arouses a windows-based error with samples of 5,000 and 250
datasets. The error is a dialogue box with the following:

"R for Windows terminal front-end has encountered a problem and needs to
close.  We are sorry for the inconvenience. If you were in the middle of
something, the information you were working on might be lost."

The new code is below. Can anyone suggest whether this error is derived
from inefficient code, or is it derived based on a windows specific
issue that can somehow be resolved and if so, how.

Thanks
Harold



library(MASS)
library(nlme)
mu<-c(100,150,200,250)
Sigma<-matrix(c(400,80,80,80,80,400,80,80,80,80,400,80,80,80,80,400),4,4
)
mu2<-c(0,0,0)
LE<-8^2 #Linking Error
Sigma2<-diag(LE,3)
sample.size<-5000
N<-100 #Number of datasets
#Take a single draw from VL distribution
vl.error<-mvrnorm(n=N, mu2, Sigma2)

intercept1 <- 0
slope1     <- 0
intercept2 <- 0
slope2     <- 0

for(i in 1:N){
temp <- data.frame(ID=seq(1:sample.size),mvrnorm(n=sample.size,
mu,Sigma))    

temp$X5 <- temp$X1
temp$X6 <- temp$X2 + vl.error[i,1] 
temp$X7 <- temp$X3 + vl.error[i,2]
temp$X8 <- temp$X4 + vl.error[i,3] 

long<-reshape(temp, idvar="ID",
varying=list(c("X1","X2","X3","X4"),c("X5","X6","X7","X8")), 
v.names=c("score.1","score.2"),direction='long')

glsrun1 <- gls(score.1~I(time-1), data=long, 
correlation=corAR1(form=~1|ID), method='ML')

glsrun2 <- gls(score.2~I(time-1), data=long, 
correlation=corAR1(form=~1|ID), method='ML')

intercept1[[i]] <- glsrun1$coefficient[1]
slope1[[i]]     <- glsrun1$coefficient[2]
intercept2[[i]] <- glsrun2$coefficient[1]
slope2[[i]]     <- glsrun2$coefficient[2]
}

cat("Sample Size =", sample.size, "\n")
cat("Number of Datasets =", N, "\n")
cat("Vertical Linking Error =", LE, "\n")
cat("Original Standard Errors","\n", "Intercept","\t",
sd(intercept1),"\n","Slope","\t","\t", sd(slope1),"\n")
cat("Modified Standard Errors","\n", "Intercept","\t",
sd(intercept2),"\n","Slope","\t","\t", sd(slope2),"\n")

rm(list=ls())
gc()

	[[alternative HTML version deleted]]


