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]]
On Thu, 20 Jan 2005 13:16:13 -0500, "Doran, Harold" <HDoran at air.org> wrote :>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.It looks to me like an nlme bug. I get the error in R-patched (built Jan 15). DrMingw shows this at the time of the crash: Rgui.exe caused an Access Violation at location 01c8ae4b in module nlme.dll Reading from location 7f1e8f18. Registers: eax=7f210020 ebx=00000000 ecx=01368c50 edx=ffffb1df esi=00004e20 edi=01108918 eip=01c8ae4b esp=0022d1d0 ebp=0022d208 iopl=0 nv up ei ng nz ac po nc cs=001b ss=0023 ds=0023 es=0023 fs=003b gs=0000 efl=00000296 Call stack: 01C8AE4B nlme.dll:01C8AE4B gls_loglik 004E5E77 R.dll:004E5E77 do_dotCode ... I changed the loop to print some status lines, and it failed after the first time it printed gls2... library(MASS) library(nlme) set.seed(123) 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){ print(i) flush.console() 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] print("reshape...") flush.console() 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') print("gls1...") flush.console() glsrun1 <- gls(score.1~I(time-1), data=long, correlation=corAR1(form=~1|ID), method='ML') print("gls2...") flush.console() 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] } Hopefully this will let someone more familiar with nlme track it down. Duncan Murdoch