[R] Failure to run mcsamp() in package arm

[Michael Kubovy]

Andrew Robinson has gently chided me for not including more  
information. So here goes:
R version 2.4.1 (2006-12-18)
powerpc-apple-darwin8.8.0

locale:
C

attached base packages:
[1] "stats"     "graphics"  "grDevices" "utils"     "methods"   "base"

other attached packages:
     foreign         car         arm   R2WinBUGS        lme4       
Matrix     lattice
    "0.8-18"     "1.2-1"    "1.0-13"     "2.0-4" "0.9975-13"  
"0.9975-11"   "0.14-16"
        MASS         JGR      iplots      JavaGD       rJava
    "7.2-32"    "1.4-15"     "1.0-5"     "0.3-5"    "0.4-14"

On Mar 7, 2007, at 4:30 PM, Michael Kubovy wrote:

> More problems. If I run
> sim(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
> from the lmer() help page.
>
> I get the error
> Error in mvrnorm(n.sims, bhat[j, ], V.beta) :
> 	'Sigma' is not positive definite
>
> On Mar 7, 2007, at 1:30 PM, Michael Kubovy wrote:
>
>> Dear r-helpers,
>>
>> I can run the examples on the mcsamp help page. For example:
>> ****************************************
>>> M1 <- lmer (y1 ~ x + (1|group))
>>> (M1.sim <- mcsamp (M1))
>> fit using lmer,
>> 3 chains, each with 1000 iterations (first 500 discarded)
>> n.sims = 1500 iterations saved
>>                            mean  sd 2.5%  25%  50%  75% 97.5% Rhat
>> n.eff
>> beta.(Intercept)           0.1 0.7 -1.2 -0.3  0.1  0.5   1.4  1.0
>> 1500
>> beta.x                     2.5 0.4  1.7  2.2  2.5  2.7   3.2  1.0
>> 1500
>> sigma.y                    3.8 0.3  3.3  3.6  3.7  3.9   4.3
>> 1.0    61
>> sigma.grop.(In)            1.5 0.8  0.0  1.0  1.4  1.9   3.3
>> 1.4    12
>> eta.group.(Intercept)[1]   0.0 1.0 -2.1 -0.5  0.0  0.6   2.0  1.0
>> 1500
>> eta.group.(Intercept)[2]   1.0 1.1 -0.9  0.2  0.9  1.7   3.4
>> 1.0    59
>> eta.group.(Intercept)[3]  -1.3 1.2 -4.0 -2.0 -1.3 -0.4   0.5
>> 1.0    66
>> eta.group.(Intercept)[4]   1.3 1.1 -0.6  0.4  1.1  2.0   3.7
>> 1.1    43
>> eta.group.(Intercept)[5]  -0.7 1.0 -3.0 -1.4 -0.6  0.0   1.2  1.0
>> 120
>> eta.group.(Intercept)[6]   1.5 1.2 -0.3  0.6  1.4  2.2   4.0
>> 1.0    49
>> eta.group.(Intercept)[7]   0.3 1.0 -1.7 -0.3  0.1  0.8   2.5  1.0
>> 440
>> eta.group.(Intercept)[8]  -1.6 1.2 -4.0 -2.4 -1.5 -0.6   0.3
>> 1.1    41
>> eta.group.(Intercept)[9]   0.4 1.0 -1.6 -0.2  0.2  0.9   2.7  1.0
>> 180
>> eta.group.(Intercept)[10] -1.0 1.1 -3.3 -1.6 -0.9 -0.2   0.8
>> 1.0    86
>>
>> For each parameter, n.eff is a crude measure of effective sample  
>> size,
>> and Rhat is the potential scale reduction factor (at convergence,
>> Rhat=1).
>> ****************************************
>> But when I try to do this with my own data I get an error:
>> ****************************************
>>> display(e7.lmer2)
>> lmer(formula = baLO ~ I(baRatio - 0.985) + delta + (1 + I(baRatio -
>> 0.985) + delta | subject), data = e7)
>>                     coef.est coef.se
>> (Intercept)        -0.19     0.06
>> I(baRatio - 0.985) -4.95     0.74
>> delta               0.41     0.06
>> Error terms:
>> Groups   Name               Std.Dev. Corr
>> subject  (Intercept)        0.13
>>            I(baRatio - 0.985) 2.57      0.45
>>            delta              0.22     -0.12 -0.94
>> Residual                    0.39
>> number of obs: 494, groups: subject, 13
>> deviance = 551.4
>>
>>> e7.sim <- mcsamp(e7.lmer2)
>> Error in as.bugs.array(sims, program = "lmer", n.iter = n.iter,
>> n.burnin = n.burnin,  :
>> 	error in parameter sigma. in parameters.to.save
>> ****************************************
>> I would appreciate a pointer to what the problem might be.
>> _____________________________
>> Professor Michael Kubovy
>> University of Virginia
>> Department of Psychology
>> USPS:     P.O.Box 400400    Charlottesville, VA 22904-4400
>> Parcels:    Room 102        Gilmer Hall
>>          McCormick Road    Charlottesville, VA 22903
>> Office:    B011    +1-434-982-4729
>> Lab:        B019    +1-434-982-4751
>> Fax:        +1-434-982-4766
>> WWW:    http://www.people.virginia.edu/~mk9y/
>>
>> ______________________________________________
>> R-help at stat.math.ethz.ch 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.
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.