[R-sig-ME] Problems with arm examples scripts
Murray Jorgensen
maj at stats.waikato.ac.nz
Thu Jan 13 04:55:58 CET 2011
I am currently running some arm scripts for some of the data sets in the
book by Andrew Gelman and Jennifer Hill. (2006). (Data Analysis Using
Regression and Multilevel/Hierarchical Models. Cambridge University Press.)
I asked a question on the R-sig-ME list and Dennis Murphy suggested
a) that I update my R to the latest version and
b) that I contact the maintainer of the arm package.
I have been mostly using a Dell laptop running Windows XP. Previously I
had been using R 2.11.1 but with packages at 2.10.1. Now I have brought
R and the packages both up to 2.12.1 but I have come across another
problem which may be an lmer problem so I am sending this email to both
R-sig-ME and Yu-Sung Su, the arm maintainer.
Here is one problem. Firstly I will give my script (which you will need
to modify to point to the location of the data file on you machine.)
library(arm)
# earnings_setup.R
he = "C:\\Files\\R&Splus\\Gelman\\examples\\earnings\\"
#library ("foreign") - arm loads this!
heights <- read.dta (paste(he,"heights.dta", sep = ""))
attach.all(heights)
age <- 90 - yearbn # survey was conducted in 1990
age[age<18] <- NA
age.category <- ifelse (age<35, 1, ifelse (age<50, 2, 3))
eth <- ifelse (race==2, 1, ifelse (hisp==1, 2, ifelse (race==1, 3, 4)))
male <- 2 - sex
ok <- !is.na (earn+height+sex+age) & earn>0 & yearbn>2
heights.clean <- as.data.frame (cbind (earn, height, sex, race, hisp,
ed, age, age.category, eth, male)[ok,])
n <- nrow (heights.clean)
# earnings_chap12.R
attach.all (heights.clean0)
# presume heights.clean intended
attach.all (heights.clean)
y <- log(earn)
x <- height
n <- length(y)
n.age <- 3
n.eth <- 4
age <- age.category
# regression of log (earnings) on height, age, and ethnicity
M1 <- lmer (y ~ x + (1 + x | eth))
display (M1)
M1 <- lmer (y ~ x + (1 + x | eth), verbose = TRUE)
(I have attached the heights.clean data as a text file to make this
example more self-contained)
In the latest version 2.12.1 lmer flags a convergence failure:
> M1 <- lmer (y ~ x + (1 + x | eth))
Warning message:
In mer_finalize(ans) : false convergence (8)
> display (M1)
lmer(formula = y ~ x + (1 + x | eth))
coef.est coef.se
(Intercept) 5.77 0.45
x 0.06 0.01
Error terms:
Groups Name Std.Dev. Corr
eth (Intercept) 0.00
x 0.00 -1.00
Residual 0.89
---
number of obs: 1187, groups: eth, 4
AIC = 3124.9, DIC = 3088.3
deviance = 3100.6
> M1 <- lmer (y ~ x + (1 + x | eth), verbose = TRUE)
0: 3113.7366: 0.0947958 0.00141417 0.00000
1: 3113.5186: 0.00000 0.00000 -0.991536
2: 3113.1254: 0.00106590 0.00111325 -0.991536
3: 3112.9224: 0.000649400 0.000651990 -0.991536
4: 3112.9133: 0.000831104 0.000668572 -0.990942
5: 3112.9121: 0.000715470 0.000730468 -0.991549
6: 3112.9117: 0.000904225 0.000536358 -0.990990
7: 3112.9114: 0.000874136 0.000514266 -0.990990
8: 3112.9112: 0.000889956 0.000518593 -0.990956
9: 3112.9111: 0.000908394 0.000472281 -0.990901
10: 3112.9110: 0.000923320 0.000460105 -0.990753
11: 3112.9106: 0.000995159 0.000309059 -0.989729
12: 3112.9103: 0.00103245 0.00000 -0.988338
13: 3112.9103: 0.00104480 1.00631e-09 -0.988338
14: 3112.9103: 0.00103952 1.72994e-08 -0.988327
15: 3112.9103: 0.00103952 2.02562e-08 -0.988327
16: 3112.9103: 0.00103951 2.01450e-08 -0.988327
17: 3112.9103: 0.00103951 2.01450e-08 -0.988327
18: 3112.9103: 0.00103951 2.01450e-08 -0.988327
Warning message:
In mer_finalize(ans) : false convergence (8)
But in my old setup 2.11.1/2.10.1 a better solution is found:
> M1 <- lmer (y ~ x + (1 + x | eth))
> display (M1)
lmer(formula = y ~ x + (1 + x | eth))
coef.est coef.se
(Intercept) 6.80 0.98
x 0.04 0.02
Error terms:
Groups Name Std.Dev. Corr
eth (Intercept) 1.37
x 0.02 -1.00
Residual 0.89
---
number of obs: 1187, groups: eth, 4
AIC = 3121.7, DIC = 3086.7
deviance = 3098.2
> M1 <- lmer (y ~ x + (1 + x | eth), verbose = TRUE)
0: 3113.7366: 0.0947958 0.00141417 0.00000
1: 3113.5186: 0.00000 0.00000 -0.991536
2: 3113.1254: 0.00106590 0.00111325 -0.991536
3: 3112.9224: 0.000649400 0.000651990 -0.991536
4: 3112.9133: 0.000831104 0.000668572 -0.990942
5: 3112.9121: 0.000715470 0.000730468 -0.991549
6: 3112.9117: 0.000904407 0.000536168 -0.990990
7: 3112.9114: 0.000874337 0.000514111 -0.990990
8: 3112.9112: 0.000890081 0.000518391 -0.990956
9: 3112.9111: 0.000908420 0.000472254 -0.990901
10: 3112.9110: 0.000923309 0.000460131 -0.990753
11: 3112.9106: 0.000994079 0.000311231 -0.989738
12: 3112.9103: 0.00103382 0.00000 -0.988215
13: 3112.9102: 0.00103963 3.55942e-10 -0.988215
14: 3112.9099: 0.00114005 0.000142779 -0.893828
15: 3112.9099: 0.00113930 9.68963e-05 -0.893828
16: 3112.9099: 0.00114352 9.70670e-05 -0.893805
17: 3112.9099: 0.00114631 9.68698e-05 -0.893782
18: 3112.9099: 0.00114593 8.96231e-05 -0.893737
19: 3112.9099: 0.00114697 8.94093e-05 -0.893635
20: 3112.9099: 0.00114840 8.74568e-05 -0.893534
21: 3112.9099: 0.00114767 8.69236e-05 -0.893331
22: 3112.9099: 0.00114960 8.45012e-05 -0.892925
23: 3112.9097: 0.00119713 0.000103484 -0.854718
24: 3112.9096: 0.00123584 0.000118803 -0.823861
25: 3112.9094: 0.00128530 0.000135095 -0.793003
26: 3112.9092: 0.00134505 0.000150925 -0.755646
27: 3112.9088: 0.00148850 0.000172432 -0.680932
28: 3112.9087: 0.00158189 0.000173348 -0.654113
29: 3112.9080: 0.00172167 0.000172246 -0.591365
30: 3112.9078: 0.00178884 0.000169040 -0.566266
31: 3112.9071: 0.00202232 2.33096e-05 -0.506894
32: 3112.9067: 0.00207199 0.00000 -0.501077
33: 3112.9063: 0.00219029 0.00000 -0.475171
34: 3112.9059: 0.00228481 0.00000 -0.456053
35: 3112.9052: 0.00250518 6.54281e-05 -0.417817
36: 3112.9048: 0.00260386 6.61956e-05 -0.402523
37: 3112.9041: 0.00281388 6.37610e-05 -0.371934
38: 3112.9037: 0.00294734 5.86692e-05 -0.358875
39: 3112.9029: 0.00315988 1.53970e-05 -0.332756
40: 3112.9023: 0.00333904 5.76185e-05 -0.315483
41: 3112.9022: 0.00337566 5.68599e-05 -0.315484
42: 3112.9022: 0.00337228 2.04001e-05 -0.315482
43: 3112.9022: 0.00339045 1.50348e-05 -0.315412
44: 3112.9022: 0.00338041 1.25433e-05 -0.315376
45: 3112.9022: 0.00338454 4.34985e-07 -0.315304
46: 3112.9007: 0.00380902 4.53446e-07 -0.277795
47: 3112.9001: 0.00401754 0.00000 -0.262792
48: 3112.8984: 0.00456608 0.00000 -0.231839
49: 3112.8956: 0.00545144 1.95070e-06 -0.195416
50: 3112.8928: 0.00652406 1.16379e-05 -0.175790
51: 3112.8916: 0.00660908 0.00000 -0.163760
52: 3112.8890: 0.00773482 4.99796e-06 -0.139726
53: 3112.8806: 0.00972168 0.00000 -0.121025
54: 3112.8786: 0.0114761 5.99051e-06 -0.0970285
55: 3112.8760: 0.0109215 6.10952e-06 -0.107449
56: 3112.8703: 0.0130487 0.00000 -0.0891341
57: 3112.8672: 0.0146986 0.00000 -0.0791535
58: 3112.8538: 0.0174766 0.00000 -0.0712682
59: 3112.8464: 0.0196425 9.24134e-08 -0.0648869
60: 3112.8229: 0.0261589 1.12351e-07 -0.0530892
61: 3112.8218: 0.0301145 0.00000 -0.0522032
62: 3112.8064: 0.0310301 2.94121e-08 -0.0482541
63: 3112.7853: 0.0370813 0.00000 -0.0428581
64: 3112.7768: 0.0392902 5.51626e-07 -0.0404835
65: 3112.7574: 0.0448577 0.00000 -0.0371558
66: 3112.7034: 0.0688045 9.08273e-06 -0.0271719
67: 3112.5859: 0.0947196 0.000137210 -0.0259377
68: 3112.5520: 0.120260 0.00000 -0.0214016
69: 3112.4267: 0.138009 0.00000 -0.0220440
70: 3112.2524: 0.241689 0.00000 -0.0175340
71: 3112.0653: 0.252030 5.74985e-06 -0.0195906
72: 3111.6754: 0.350688 0.00000 -0.0176928
73: 3111.5729: 0.396312 0.00000 -0.0169346
74: 3111.4491: 0.441941 0.00000 -0.0166439
75: 3110.9290: 0.624460 0.00000 -0.0161161
76: 3110.0355: 1.03912 0.00000 -0.0160134
77: 3109.9542: 1.12017 0.00000 -0.0157653
78: 3109.8078: 1.29685 0.00000 -0.0157324
79: 3109.7519: 1.45344 0.00000 -0.0157254
80: 3109.7469: 1.52828 0.00000 -0.0156867
81: 3109.7461: 1.53677 0.00000 -0.0157060
82: 3109.7460: 1.53350 0.00000 -0.0157010
83: 3109.7460: 1.53352 0.00000 -0.0157010
The two iterations start off at the same values but start to diverge
around step 12.
Murray
--
Dr Murray Jorgensen http://www.stats.waikato.ac.nz/Staff/maj.html
Department of Statistics, University of Waikato, Hamilton, New Zealand
Email: maj at waikato.ac.nz Fax 7 838 4155
Phone +64 7 838 4773 wk Home +64 7 825 0441 Mobile 021 0200 8350
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: heights.clean.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110113/413ab686/attachment-0004.txt>
More information about the R-sig-mixed-models
mailing list