[R] NLME questions -- interpretation of results

ctu at bigred.unl.edu ctu at bigred.unl.edu
Thu Jul 3 06:44:25 CEST 2008


Hi Jenny, (I use the data you provide in the previous e-mail)
For the 1st question, let me assume you only want to compare loc: A vs. B
So you could specified your code like this:
fmAB <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST,
             random = Asym ~1,
             fixed = Asym+R0+lrc ~ loc,
             start=c(0.97,0,
                     1.14,0,
                    -0.18,0))
> summary(fmAB)
Nonlinear mixed-effects model fit by maximum likelihood
   Model: Y ~ SSlogis(X, Asym, R0, lrc)
  Data: LAST
         AIC       BIC   logLik
   -31.02303 -19.70017 24.51152

Random effects:
  Formula: Asym ~ 1 | loc
         Asym.(Intercept)
StdDev:     1.549779e-06

  Formula: Asym ~ 1 | dir %in% loc
         Asym.(Intercept)   Residual
StdDev:     6.124237e-08 0.09426086

Fixed effects: Asym + R0 + lrc ~ loc
                       Value Std.Error DF   t-value p-value
Asym.(Intercept)  0.9477404 0.1037337 14  9.136280  0.0000
Asym.locB         0.0982456 0.4175971 14  0.235264  0.8174
R0.(Intercept)    1.1289387 0.0652189 14 17.310003  0.0000
R0.locB           0.1390656 0.3946717 14  0.352358  0.7298
lrc.(Intercept)  -0.2110057 0.0656513 14 -3.214036  0.0062
lrc.locB         -0.0820484 0.6826382 14 -0.120193  0.9060

Then you know Asym, R0, and lrc of loc B are not significant.  
Moreover, you can test the joint fixed effect by
anova(fmAB)(Pinherio and Bate, 2000 Book, p 374)

for the 2nd question, How to get the fitted value for particular level?
Based on this example, let me assume you want to get the fitted value of A/N.

then you could write a small code like this:
> FV<-data.frame(F.V=fitted(fmAB), group=summary(fmAB)$groups$dir)
> A.N<-FV[is.element(FV$group, c("A/N")),]
          F.V group
9  0.4209011   A/N
10 0.2726129   A/N

hope this is helpful~

Chunhao












Quoting Jenny Sun <jenny.sun.sun at gmail.com>:

> Thank you for your reply Chunhao!
>
> I attached only part of the test data and that is why you might not   
> be able to get convergence. Sorry.
>
> I  have a couple more questions:
>
> For the second question you answered, how to specify the correct   
> length of starting values. I tried using the length of levels in   
> each of the parameters in the start list but found:
>
>> fm1 <- nlme(DIFN ~ SSlogis(SVA, Asym, R0, lrc),data = LAST,fixed =   
>> Asym + R0 + lrc ~ dir %in% loc,random = Asym ~ 1,start =list(Asym =  
>>  c(1,1,1,1), R0 =  c(1,1,1,1), lrc =  c(-5,-2,-2,-2)))
> Error in nlme.formula(DIFN ~ SSlogis(SVA, Asym, R0, lrc), data = LAST,  :
>   start must have a component called "fixed"
>
> I've got two loc levels (A,B) with four group levels(N,E,S,W);   How  
>  I am gonna define the list and the component called"fixed"?
>
> My another question is about the fitted value of the model. If I   
> want to calculate adjusted R square, I have to get fitted(fm1).   
> WHich has values like this;
>
>  >fitted(fm1)[1:40]
>      AB/N      AB/N      AB/E      AB/S      AB/W      AB/W        
> AB/W      AB/W      AB/W      AB/W
> 0.6541876 0.7421748 0.8408251 0.5879220 0.4889387 0.6129576   
> 0.5097593 0.6195679 0.5152567 0.5680860
>      AB/W      AB/W      AB/W      AB/W      AB/W      AB/W        
> AB/N      AB/N      AB/N      AB/E
> 0.4724423 0.8128148 0.7674529 0.7106698 0.6553155 0.6074771   
> 0.5036201 0.5464105 0.6062978 0.6878438
>      AB/N      AB/N      AB/N      AB/S      AB/S      AB/S        
> AB/S      AB/S      AB/S      AB/S
> 0.7792725 0.8411961 0.7942503 0.7354845 0.5895700 0.6781973   
> 0.6286886 0.5212052 0.8864748 0.8370021
>      AB/S      AB/S      AB/N      AB/N      AB/N      AB/N        
> AB/E      AB/E      AB/E      AB/E
> 0.7750731 0.7147024 0.6625288 0.5492599 0.5959280 0.6612426   
> 0.7501786 0.8498928 0.6274681 0.7118615
>
> My question is how to get the fitted values for specified group   
> levels (eg. values for AB/E)?
>
> Thank all very much!
>
> Jenny
>
>
>> Hi Jenny,
>> I try your code but I did not get in converge in fm3 (see the below).
>> For the first question, you could use fm1 to interpret the result
>> without bothering fm2 and fm3. It means that R0 and lrc can be treated
>> as pure fixed effects (Pinherir and Bates, 2000 Book).
>>
>> 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~
>>
>> Hope this helpful.
>>
>> Chunhao Tu
>>
>>> test<-read.table(file="C:\\Documents and
>>> Settings\\ado_cabgfaculty\\Desktop\\sun.txt", header=T)
>>> 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
>>>
>>>
>
>



More information about the R-help mailing list