Yu-Kang -
Simulations by their nature use randomly generated data.
Sometimes the random data doesn't contain enough information
to fully determine the parameter estimates for one iteration
or another. It seems likely that that is what happened here.
The design matrix is singular for one iteration (maybe there
are NO simulated subjects in one arm of the trial) and
backsolve() very properly returns an error message. On the
very next iteration, with a different randomly-generated data
set, everything should work, again.
The function try() allows you to get past these problematic
iterations. It's better, of course, to figure out what the
requirements for a fully-specified data set are, and arrange
the random generator so that these are guaranteed to be met.
- tom blackwell - u michigan medical school - ann arbor -
On Tue, 28 Oct 2003, Tu Yu-Kang wrote:
> Dear R-users,
>
> I am a dentist (so forgive me if my question looks stupid) and came across
> a problem when I did simulations to compare a few single level and two
> level regressions.
>
> The simulations were interrupted and an error message came out like
'Error
> in MEestimate(lmeSt, grps) : Singularity in backsolve at level 0, block
1'.
>
> My collegue suggested that this might be due to my codes was not efficient
> and ran out of memory. If this is the reason, could you please help me
> improve my codes writing.
>
> However, as I slightly changed the parameters, it ran well. So I suspect
> memory is problem.
>
> I use R 1.8.0 and Windows XP professional. My computer has a Pentium 4 2.4
> with 512 MB memory.
>
> Thanks in advance.
>
> best regards,
>
> Yu-Kang Tu
>
> Clinical Research Fellow
> Leeds Dental Institute
> University of Leeds
>
> ## change scores simulation
> close.screen(all=TRUE)
> split.screen(c(3,3))
> nitns<-10000
> nsims<-100
> r<-0.1
> param1<-c(1:nitns)
> param2<-c(1:nitns)
> param3<-c(1:nitns)
> param4<-c(1:nitns)
> param5<-c(1:nitns)
> param6<-c(1:nitns)
> param7<-c(1:nitns)
> param8<-c(1:nitns)
> param9<-c(1:nitns)
> param10<-c(1:nitns)
> param11<-c(1:nitns)
> param12<-c(1:nitns)
> param13<-c(1:nitns)
> param14<-c(1:nitns)
> for(itn in 1:nitns){
> g<-rbinom(nsims,1,0.5)
> b<-rnorm(nsims,0,1)*10
> rn<-rnorm(nsims,0,1)*10
> a<-b*r+rn*(1-r^2)^0.5
> a<-round(a)+50
> a<-a-g*5
> b<-round(b)+50
> abs.2<-function(x) ifelse(x<1,1,x)
> b<-abs.2(b)
> c<-b-a
> p<-c/b
> lm1<-lm(a~g)
> lm2<-lm(c~g)
> lm3<-lm(p~g)
> lm4<-lm(a~b+g)
> gr<-c(g,g)
> occasion<-rep(0:1,c(nsims,nsims))
> occ<-occasion-0.5
> ppd<-c(b,a)
> h<-rep(0,nsims)
> mb<-mean(b)
> bppd<-b-mb
> bappd<-c(h,bppd)
> occgr<-occ*gr
> subject<-c(1:nsims)
> sub<-c(subject,subject)
> library(nlme)
> lm5<-lme(ppd~occ+gr+occgr,random=~1|sub)
> lm5f<-fixed.effects(lm5)
> lm5c<-as.matrix(lm5f)
> lm5a<-anova(lm5)
> lm6<-lme(ppd~occ+gr+occgr,random=~1|sub,method="ML")
> lm6f<-fixed.effects(lm6)
> lm6c<-as.matrix(lm6f)
> lm6a<-anova(lm6)
> lm7<-lme(ppd~occ+gr+occgr+bappd,random=~1|sub,method="ML")
> lm7f<-fixed.effects(lm7)
> lm7c<-as.matrix(lm7f)
> lm7a<-anova(lm7)
> param1[itn]<-coef(summary(lm1))[2,1]
> param2[itn]<-coef(summary(lm2))[2,1]
> param3[itn]<-coef(summary(lm3))[2,1]
> param4[itn]<-coef(summary(lm4))[3,1]
> param5[itn]<-lm5c[4,1]
> param6[itn]<-lm6c[4,1]
> param7[itn]<-lm7c[4,1]
> param8[itn]<-coef(summary(lm1))[2,4]
> param9[itn]<-coef(summary(lm2))[2,4]
> param10[itn]<-coef(summary(lm3))[2,4]
> param11[itn]<-coef(summary(lm4))[3,4]
> param12[itn]<-lm5a[4,4]
> param13[itn]<-lm6a[4,4]
> param14[itn]<-lm7a[4,4]
> }
> #the error message came out here.
> #But if I change some of the variables to:
> b<-rnorm(nsims,0,1)*2
> rn<-rnorm(nsims,0,1)*2
> a<-b*r+rn*(1-r^2)^0.5
> a<-round(a)+7
> a<-a-g*2
> b<-round(b)+9
> abs.1<-function(x) ifelse(x<5,5,x)
> b<-abs.1(b)
> abs.2<-function(x) ifelse(x<1,1,x)
> a<-abs.2(a)
>
> There is no poblem to complete the simulations.
>
> _________________________________________________________________
> ²{¦b´N¤W MSN ¥xÆWºô¯¸¡G»P¿ËªB¦n¤Íºò±KÁpô¡A§Y®É´x´¤·s»D¡B°]¸g¡B®T¼Öªº³Ì·s°T
> ®§ http://msn.com.tw
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>