R help - Oct 2009 - Interpreting gnls() output in comparison to nls()

Hi,

I've been trying to work with the gnls() function in the "nlme"
package.  My
decision to use gnls() was so that I could fit varPower and such to some of 
the data.  However, in working with a small dataset, I've found that the 
results given by gnls() don't seem to make any sense and they differ 
substantially from those produced by nls().  I suspect that I am just 
misinterpreting the gnls() output and could use a little guidance.


Here's the data I am trying to fit:
------------------------------------
> myPbmcData
Grouped Data: lnCount ~ Time | Type
    Time Mouse Count   Type  lnCount
1     0  T0-1  37.6  Naive 3.627004
2     0  T0-2  23.6  Naive 3.161247
3     0  T0-3  49.2  Naive 3.895894
4     8  T8-1  20.8  Naive 3.034953
5     8  T8-2  26.9  Naive 3.292126
6     8  T8-3  34.0  Naive 3.526361
7    24 T24-1  36.8  Naive 3.605498
8    24 T24-2  34.0  Naive 3.526361
9    24 T24-3  19.6  Naive 2.975530
10   48 T48-1  55.4  Naive 4.014580
11   48 T48-2  54.2  Naive 3.992681
12   48 T48-3  51.4  Naive 3.939638
13    0  T0-1 342.0 Memory 5.834811
14    0  T0-2 139.0 Memory 4.934474
15    0  T0-3 191.0 Memory 5.252273
16    8  T8-1 104.0 Memory 4.644391
17    8  T8-2 192.0 Memory 5.257495
18    8  T8-3 204.0 Memory 5.318120
19   24 T24-1 161.0 Memory 5.081404
20   24 T24-2 131.0 Memory 4.875197
21   24 T24-3 230.0 Memory 5.438079
22   48 T48-1 458.0 Memory 6.126869
23   48 T48-2 594.0 Memory 6.386879
24   48 T48-3 374.0 Memory 5.924256
>
------------------------------------
Now I fit nls() and gnls() to the data.

--------------------------------------

##Fit nls
> nlsOut +   nls(lnCount ~ log(C0[Type] + C1[Type] * Time + C2[Type]*
Time^2),
+       start = list( C0 = c(exp(3.5), exp(3.5)), C1 = c(-0.01, -0.01), C2 = 
c(0.01, 0.01)),
+        data = myPbmcData,
+       )
>
> summary(nlsOut)
Formula: lnCount ~ log(C0[Type] + C1[Type] * Time + C2[Type] * Time^2)

Parameters:
      Estimate Std. Error t value Pr(>|t|)
C01  33.82469    5.08870   6.647 3.08e-06 ***
C02 209.40868   31.01673   6.751 2.51e-06 ***
C11  -0.90447    0.58030  -1.559  0.13649
C12  -8.64779    3.73034  -2.318  0.03241 *
C21   0.02775    0.01261   2.201  0.04102 *
C22   0.29156    0.08975   3.249  0.00446 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1

Residual standard error: 0.3 on 18 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 1.447e-06
>
>
> ##Fit gnls
> gnlsOut +   gnls(lnCount ~ log(C0 + C1 * Time + C2* Time^2),
+       start = list( C0 = c(exp(3.5), exp(3.5)), C1 = c(-0.5, -0.5), C2 = 
c(0.01, 0.01)),
+        data = myPbmcData,
+        params = list( C0 + C1 + C2 ~ Type),
+        verbose = TRUE
+     )
>
> summary(gnlsOut)
Generalized nonlinear least squares fit
   Model: lnCount ~ log(C0 + C1 * Time + C2 * Time^2)
   Data: myPbmcData
        AIC      BIC    logLik
   17.41510 25.66148 -1.707552

Coefficients:
                    Value Std.Error   t-value p-value
C0.(Intercept) 121.61674 15.715703  7.738549  0.0000
C0.Type.L      124.15665 22.225361  5.586260  0.0000
C1.(Intercept)  -4.77613  1.887602 -2.530265  0.0209
C1.Type.L       -5.47536  2.669472 -2.051104  0.0551
C2.(Intercept)   0.15965  0.045315  3.523169  0.0024
C2.Type.L        0.18654  0.064085  2.910877  0.0093

  Correlation:
                C0.(I) C0.T.L C1.(I) C1.T.L C2.(I)
C0.Type.L       0.948
C1.(Intercept) -0.750 -0.713
C1.Type.L      -0.713 -0.750  0.953
C2.(Intercept)  0.554  0.529 -0.924 -0.884
C2.Type.L       0.529  0.554 -0.884 -0.924  0.961

Standardized residuals:
         Min          Q1         Med          Q3         Max
-1.41262246 -0.76633639 -0.03330171  0.67865673  1.63503742

Residual standard error: 0.300007
Degrees of freedom: 24 total; 18 residual
>
>
--------------------------------------------
Examining the the results, they don't match up as I read them.
For example, look at C0. nls() gives initial values for
     CO Naive ~33 and  CO Memory ~210.

In contrast, gnls() gives an intercept at C0 121 and an effect of, if I am 
reading the output correctly,
    C0 Naive ~ 121.62- 124.16 = -2.54 and CO Memory ~ 121.62+ 124.16 = 245.78

However, if I compare the predictions they match up very well

--------------------------------------------------
> (predict(gnlsOut) - predict(nlsOut))/predict(nlsOut)
             1             2             3             4             5
  3.646201e-07  3.646201e-07  3.646201e-07  7.849412e-07  7.849412e-07
             6             7             8             9            10
  7.849412e-07  3.655158e-07  3.655158e-07  3.655158e-07 -1.298393e-06
            11            12            13            14            15
-1.298393e-06 -1.298393e-06  6.636562e-08  6.636562e-08  6.636562e-08
            16            17            18            19            20
-7.458066e-10 -7.458066e-10 -7.458066e-10 -7.685896e-08 -7.685896e-08
            21            22            23            24
-7.685896e-08  1.456831e-08  1.456831e-08  1.456831e-08
attr(,"label")
[1] "Fitted values"
---------------------------------------------------
Could someone set me straight as to how I should be interpreting the gnls() 
output?


Many thanks for your time.

Mike


-----------------------------------------------------
Department of Ecology & Evolutionary Biology
569 Dabney Hall
University of Tennessee
Knoxville, TN 37996-1610

phone:(865) 974-6453
fax:  (865) 974-6042

web: http://eeb.bio.utk.edu/gilchrist.asp