> library(nlme) > help(Soybean) #Soybean leaf weight over time, described in Davidian and Giltinan (1995, 1.1.3, p.7) #from an experiment to compare growth patterns of two genotypes of soybeans: #Plant Introduction #416937 (P), an experimental strain, and Forrest (F), a commercial variety. > Soybean[1:3,] Grouped Data: weight ~ Time | Plot Plot Variety Year Time weight 1 1988F1 F 1988 14 0.106 2 1988F1 F 1988 21 0.261 3 1988F1 F 1988 28 0.666 > plot(Soybean) # average leaf weight per plant of 2 types, over 3 planting years > plot(Soybean, outer=~Year*Variety) #fit nonlinear model > fit <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), data=Soybean) > summary(fit) Formula: weight ~ SSlogis(Time, Asym, xmid, scal) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 18.4182 0.4363 42.22 <2e-16 *** xmid 53.9562 0.6900 78.20 <2e-16 *** scal 8.1535 0.4671 17.45 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.264 on 409 degrees of freedom #could use nlme with starting values from nls fit and following help file for nlme #this is the same as fm1.nlme (below) fm1b.nlme <- nlme(weight~Asym/(1+exp((xmid-Time)/scal)), data=Soybean, fixed = Asym + xmid + scal ~ 1, random = Asym + xmid + scal ~ 1, start=c(Asym=18,xmid=53,scal=8)) #A nicer way is to use: #nlsList function fits nonlinear model to each profile > summary(fm1 <- nlsList(weight ~ SSlogis(Time, Asym, xmid, scal), data = Soybean)) Error in qr.solve(QR.B, cc) : singular matrix 'a' in solve Call: Model: weight ~ SSlogis(Time, Asym, xmid, scal) | Plot Data: Soybean Coefficients: Asym Estimate Std. Error t value Pr(>|t|) 1988F4 15.151570 0.8934694 16.9581296 4.463395e-07 1988F2 19.745419 1.9870782 9.9369110 3.794038e-07 ... 1989P8 NA NA NA NA ... 1990P2 11.656903 2.101629 5.546605 2.695014e-05 1990P4 10.971568 1.844637 5.947820 4.695505e-03 Residual standard error: 1.043835 on 263 degrees of freedom #Error message from nls indicates convergence was not attained for #plot 1989P8 #Now, fit random effects for all 3 parameters in the logistic (Asym, xmid, scal) #in the nonlinear mixed effects model, fixed effects are the mean values of the #parameters and random effects are deviations from their mean values. > fm1.nlme <- nlme(fm1) > summary(fm1.nlme) Nonlinear mixed-effects model fit by maximum likelihood Model: weight ~ SSlogis(Time, Asym, xmid, scal) Data: Soybean AIC BIC logLik 1499.671 1539.881 -739.8354 Random effects: Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Level: Plot Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr Asym 5.201174 Asym xmid xmid 4.197457 0.721 scal 1.404713 0.711 0.958 Residual 1.123463 Fixed effects: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Value Std.Error DF t-value p-value Asym 19.25302 0.8031906 362 23.97067 0 xmid 55.01985 0.7272711 362 75.65246 0 scal 8.40332 0.3152871 362 26.65292 0 Correlation: Asym xmid xmid 0.724 scal 0.620 0.807 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -6.08689827 -0.22161218 -0.03390493 0.29740625 4.84685692 Number of Observations: 412 Number of Groups: 48 > plot(fm1.nlme) #residuals show increasing variance for within group errors (heteroscedasticity) #Can use the power variance function, varPower, Var(e_ij) = sigma^2*|v_ij|^{2*delta} #see help(varPower) > fm2.nlme <- update(fm1.nlme, weights=varPower()) > anova(fm1.nlme, fm2.nlme) Model df AIC BIC logLik Test L.Ratio p-value fm1.nlme 1 10 1499.6709 1539.8811 -739.8354 fm2.nlme 2 11 745.7111 789.9423 -361.8555 1 vs 2 755.9598 <.0001 > anova(fm1.nlme, fm2.nlme) > summary(fm2.nlme) > summary(fm2.nlme) Nonlinear mixed-effects model fit by maximum likelihood Model: weight ~ SSlogis(Time, Asym, xmid, scal) Data: Soybean AIC BIC logLik 745.7111 789.9423 -361.8555 Random effects: Formula: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Level: Plot Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr Asym 2.8923715276 Asym xmid xmid 1.2871418773 -1.000 scal 0.0000789892 0.021 -0.021 Residual 0.2267566858 Variance function: Structure: Power of variance covariate Formula: ~fitted(.) Parameter estimates: power 0.9491257 Fixed effects: list(Asym ~ 1, xmid ~ 1, scal ~ 1) Value Std.Error DF t-value p-value Asym 16.93426 0.5577296 362 30.36285 0 xmid 51.81873 0.4424837 362 117.10880 0 scal 7.53841 0.0859232 362 87.73430 0 Correlation: Asym xmid xmid 0.155 scal 0.310 0.769 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.5816748 -0.6134148 -0.1105220 0.5885103 3.4706660 Number of Observations: 412 Number of Groups: 48 #better fit! > plot(fm2.nlme) > qqnorm(fm2.nlme)