> dat.long = read.csv("http://www.ics.uci.edu/~staceyah/112-203/data/autism.csv") > head(dat.long) age vsae sicdegp childid 1 2 6 3 1 2 3 7 3 1 3 5 18 3 1 4 9 25 3 1 5 13 27 3 1 6 2 17 3 3 > dat.long$sicdegp = factor(dat.long$sicdegp) > dat.long$childid = factor(dat.long$childid) > library(lattice) > xyplot(vsae ~ age | childid, data=dat.long, type="l") > library(nlme) > ?lme > mod1 = lme(vsae ~ sicdegp*age, random = ~ 1 | childid, data = dat.long, + weights = varIdent(form = ~1|age)) Error in na.fail.default(list(age = c(2L, 3L, 5L, 9L, 13L, 2L, 3L, 5L, : missing values in object > > unique(dat.long$sicdegp) [1] 3 2 1 Levels: 1 2 3 > mod1 = lme(vsae ~ sicdegp*age, random = ~ 1 | childid, data = dat.long, na.action=na.omit, + weights = varIdent(form = ~1|age)) > summary(mod1) Linear mixed-effects model fit by REML Data: dat.long AIC BIC logLik 4629.539 4682.382 -2302.769 Random effects: Formula: ~1 | childid (Intercept) Residual StdDev: 2.900144 1.642413 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | age Parameter estimates: 2 3 5 9 13 1.000000 4.011672 6.267865 16.883237 26.783182 Fixed effects: vsae ~ sicdegp * age Value Std.Error DF t-value p-value (Intercept) 1.141038 0.9138404 449 1.248618 0.2125 sicdegp2 0.055366 1.2092679 155 0.045784 0.9635 sicdegp3 -3.448228 1.3387987 155 -2.575613 0.0109 age 2.973609 0.3705607 449 8.024621 0.0000 sicdegp2:age 0.822015 0.4898773 449 1.678001 0.0940 sicdegp3:age 4.405273 0.5364819 449 8.211411 0.0000 Correlation: (Intr) scdgp2 scdgp3 age scdg2: sicdegp2 -0.756 sicdegp3 -0.683 0.516 age -0.859 0.649 0.587 sicdegp2:age 0.650 -0.859 -0.444 -0.756 sicdegp3:age 0.594 -0.449 -0.853 -0.691 0.522 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.4417142 -0.5616347 -0.1070310 0.3404212 5.9675233 Number of Observations: 610 Number of Groups: 158 > with(dat.long, interaction.plot(age,childid,vsae,fun=function(x){mean(x,na.rm=T)},trace.label=FALSE)) > mod2 = lme(vsae ~ sicdegp*age, random = ~ age | childid, data = dat.long, na.action=na.omit) Error in lme.formula(vsae ~ sicdegp * age, random = ~age | childid, data = dat.long, : nlminb problem, convergence error code = 1 message = iteration limit reached without convergence (10) > ?lme > mod2 = lme(vsae ~ sicdegp*age, random = ~ age | childid, data = dat.long, na.action=na.omit, + weights = varIdent(form = ~1|age)) > anova(mod1,mod2) Model df AIC BIC logLik Test L.Ratio mod1 1 12 4629.539 4682.382 -2302.769 mod2 2 14 4479.361 4541.011 -2225.680 1 vs 2 154.178 p-value mod1 mod2 <.0001 > mod3 = lme(vsae ~ sicdegp*age, random = ~ age | childid, data = dat.long, na.action=na.omit, + weights = varIdent(form = ~1)) Error in lme.formula(vsae ~ sicdegp * age, random = ~age | childid, data = dat.long, : nlminb problem, convergence error code = 1 message = iteration limit reached without convergence (10) > ?lme > summary(mod2) Linear mixed-effects model fit by REML Data: dat.long AIC BIC logLik 4479.361 4541.011 -2225.68 Random effects: Formula: ~age | childid Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 6.563872 (Intr) age 3.476953 -0.948 Residual 2.422935 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | age Parameter estimates: 2 3 5 9 13 1.000000 2.312368 3.346606 5.478704 6.863456 Fixed effects: vsae ~ sicdegp * age Value Std.Error DF t-value p-value (Intercept) 1.117675 1.1515112 449 0.970616 0.3323 sicdegp2 0.357681 1.5221662 155 0.234982 0.8145 sicdegp3 -3.019978 1.6923876 155 -1.784448 0.0763 age 3.079967 0.5545747 449 5.553745 0.0000 sicdegp2:age 0.624896 0.7327183 449 0.852846 0.3942 sicdegp3:age 4.180156 0.8151542 449 5.128056 0.0000 Correlation: (Intr) scdgp2 scdgp3 age scdg2: sicdegp2 -0.756 sicdegp3 -0.680 0.515 age -0.922 0.697 0.627 sicdegp2:age 0.698 -0.922 -0.475 -0.757 sicdegp3:age 0.627 -0.474 -0.919 -0.680 0.515 Standardized Within-Group Residuals: Min Q1 Med Q3 -2.58428338 -0.42192951 -0.07454467 0.39123757 Max 6.09939992 Number of Observations: 610 Number of Groups: 158 > mod2 = lme(vsae ~ sicdegp*age, random = ~ age | childid, data = dat.long, na.action=na.omit, + weights = varIdent(form = ~1|age)) > summary(mod2) Linear mixed-effects model fit by REML Data: dat.long AIC BIC logLik 4479.361 4541.011 -2225.68 Random effects: Formula: ~age | childid Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 6.563872 (Intr) age 3.476953 -0.948 Residual 2.422935 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | age Parameter estimates: 2 3 5 9 13 1.000000 2.312368 3.346606 5.478704 6.863456 Fixed effects: vsae ~ sicdegp * age Value Std.Error DF t-value p-value (Intercept) 1.117675 1.1515112 449 0.970616 0.3323 sicdegp2 0.357681 1.5221662 155 0.234982 0.8145 sicdegp3 -3.019978 1.6923876 155 -1.784448 0.0763 age 3.079967 0.5545747 449 5.553745 0.0000 sicdegp2:age 0.624896 0.7327183 449 0.852846 0.3942 sicdegp3:age 4.180156 0.8151542 449 5.128056 0.0000 Correlation: (Intr) scdgp2 scdgp3 age scdg2: sicdegp2 -0.756 sicdegp3 -0.680 0.515 age -0.922 0.697 0.627 sicdegp2:age 0.698 -0.922 -0.475 -0.757 sicdegp3:age 0.627 -0.474 -0.919 -0.680 0.515 Standardized Within-Group Residuals: Min Q1 Med Q3 -2.58428338 -0.42192951 -0.07454467 0.39123757 Max 6.09939992 Number of Observations: 610 Number of Groups: 158 > mod3 = lme(vsae ~ sicdegp*age, random = ~ age | childid, data = dat.long, na.action=na.omit, + weights = varIdent(form = ~1)) Error in lme.formula(vsae ~ sicdegp * age, random = ~age | childid, data = dat.long, : nlminb problem, convergence error code = 1 message = iteration limit reached without convergence (10)