# Re: [R] NLME questions -- interpretation of results

From: <ctu_at_bigred.unl.edu>
Date: Wed, 02 Jul 2008 17:08:22 -0500

For the second question, your want to know "is AB/E different from the AB/S"

The simplest way is to change your fixed statement: fixed = Asym+R0+lrc ~ dir %in% loc
and specify the correct length of starting values.

If I am wrong please correct me~

Chunhao Tu

> LAST<-groupedData(Y~X|loc/dir, data=test)
>
> fm1 <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST,

```+             random = Asym ~1,
+             fixed = Asym+R0+lrc ~ 1,
+             start=c(Asym = 0.97, R0 =  1.14, lrc =  -0.18))

> fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1))

```
> fm3 <- update(fm2, random = pdDiag(Asym+R0+lrc~ 1))
Error in nlme.formula(model = Y ~ SSlogis(X, Asym, R0, lrc), data = LAST, :

Step halving factor reduced below minimum in PNLS step

Quoting Jenny Sun <jenny.sun.sun_at_gmail.com>:

> My special thanks to Chunhao Tu for the suggestions about testing
> significance of two locations.
>
> I used logistic models to describe relationships between Y and X at
> two locations (A & B). And within each location, I have four groups
> (N,E,S,W)representing directions. So the test data can be arranged as:
>
> Y X dir loc
> 0.6295 0.8667596 S A
> 0.7890 0.7324820 S A
> 0.4735 0.9688875 S A
> 0.7805 1.1125239 S A
> 0.8640 0.9506174 E A
> 0.9445 0.6582157 E A
> 0.8455 0.5558860 E A
> 0.9380 0.3304870 E A
> 0.4010 1.1763090 N A
> 0.2585 1.3202890 N A
> 0.3750 1.1763090 E A
> 0.3855 1.3202890 E A
> 0.3020 1.1763090 S A
> 0.2300 1.3202890 S A
> 0.3155 1.1763090 W A
> 0.8890 0.6915861 W B
> 0.9185 0.6149019 W B
> 0.9275 0.5289258 W B
> 0.8365 0.9507088 S B
> 0.7720 0.8842165 N B
> 0.8615 0.8245123 N B
> 0.9170 0.7559687 W B
> 0.9590 0.6772720 W B
> 0.9900 0.5872023 W B
> 0.9940 0.4849064 W B
> 0.7500 0.9560776 W B
>
>
> The data is grouped using:
>
>> LAST<-groupedData(Y~X|loc/dir, data=test)
>
> I then used logistic models to define the relationship between Y and
> X, and got fm1, fm2, and fm3 as follows:
>
> --------------------------
>> fm1 <- nlme(DIFN ~ SSlogis(SVA, Asym, R0, lrc),data = LAST,fixed =
>> Asym + R0 + lrc ~ 1,random = Asym ~ 1,start =c(Asym = 1, R0 = 1,
>> lrc = -5))
>> fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1))
>> fm3 <- update(fm2, random = pdDiag(Asym + R0 + lrc ~ 1))
>> anova(fm1,fm2,fm3)
> ------------------------------------------------------------
> ANOVA showed:
>
>> anova(fm1,fm2,fm3)
> Model df AIC BIC logLik Test L.Ratio p-value
> fm1 1 7 -1809.913 -1774.304 910.9564
> fm2 2 9 -1805.774 -1758.295 910.8871 1 vs 2 0.1386696 0.9999
> fm3 3 12 -1801.822 -1742.473 910.9109 2 vs 3 0.0475543 0.9666
>
> ** question: do the results show that fm1 could represent the
> results of fm2 and fm3?
>
>> coef(fm1)
> Asym R0 lrc
> AB/E 0.9148927 1.389432 -0.3009858
> AB/N 0.8775250 1.389432 -0.3009858
> AB/S 0.9247592 1.389432 -0.3009858
> AB/W 0.8479180 1.389432 -0.3009858
> BC/E 0.8791908 1.389432 -0.3009858
> BC/N 0.8414229 1.389432 -0.3009858
> BC/S 0.9169323 1.389432 -0.3009858
> BC/W 0.8817838 1.389432 -0.3009858
>
> ** question: how could I know if any of the models is significantly
> different from the other ones? (eg. AB/E is different from the AB/S)?
>
>> summary(fm1)
> Nonlinear mixed-effects model fit by maximum likelihood
> Model: DIFN ~ SSlogis(SVA, Asym, R0, lrc)
> Data: LAST
> AIC BIC logLik
> -1809.913 -1774.304 910.9564
>
> Random effects:
> Formula: Asym ~ 1 | loc
> Asym
> StdDev: 2.303402e-05
>
> Formula: Asym ~ 1 | dir %in% loc
> Asym Residual
> StdDev: 0.03208693 0.1741559
>
> Fixed effects: Asym + R0 + lrc ~ 1
> Value Std.Error DF t-value p-value
> Asym 0.8855531 0.015375906 2783 57.59355 0
> R0 1.3894322 0.009418047 2783 147.52869 0
> lrc -0.3009858 0.012833066 2783 -23.45393 0
> Correlation:
> Asym R0
> R0 -0.440
> lrc -0.452 0.150
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -4.1326757 -0.6117037 0.1082112 0.6575250 3.3297270
>
> Number of Observations: 2793
> Number of Groups:
> loc dir %in% loc
> 2 8
>
>
> I have marked all the codes and questions(**). Any answers and
> suggestions are appreciated.
>
> Have a good day!
>
> Jenny
>
>

R-help_at_r-project.org 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. Received on Wed 02 Jul 2008 - 22:11:34 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 03 Jul 2008 - 02:31:44 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.