R help - May 2008 - How to test significant differences for non-linear relationships for two locations

Hi List,

I have to compare a relationship between y and x for two locations. I found
logistic regression fits both datasets well, but I am not sure  how to test if
relationships for both sites are significantly different. I searched the r site,
however no answers exactly match the question.

I used Tukey's HSD to compare two means, but the relationship in my study
was not simply linear. So I was wondering if there is anyone had experience in
making such comparisons. Thanks in advance!

Jenny

	[[alternative HTML version deleted]]

Hi Jenny,

My suggestion is that
1. require(nlme)
2. combine two data sets and create a new variable e.g., location
thus your new data set will look like this as following, such as,

y    x location
10   0   L1
11   1   L1
12   2   L1
10.1 0   L2
11.2 1   L2
12.1 2   L3

3. Group your data set, LAST<-groupedData(y~x|location, data=new)
4. test<-nlme(y~logistic.function, data=LAST, fixed=alpha+tau+psi~location,
               random=pdDiag(alpha+tau+psi~1),
               start=c(alpha.value,0, tau.value,0,psi.valu,0))
5.anova(test) (<-your answer)
You could find this in Pinheiro and Bates books

Chunhao Tu


Quoting Jenny Sun <jenny.sun.sun at gmail.com>:
> Hi List,
>
> I have to compare a relationship between y and x for two locations.   
> I found logistic regression fits both datasets well, but I am not   
> sure  how to test if relationships for both sites are significantly   
> different. I searched the r site, however no answers exactly match   
> the question.
>
> I used Tukey's HSD to compare two means, but the relationship in my   
> study was not simply linear. So I was wondering if there is anyone   
> had experience in making such comparisons. Thanks in advance!
>
> Jenny
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> 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.
>

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

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
>
>