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
______________________________________________
R-help@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.
