[R] Grouped Regression

Bill.Venables@CMIS.CSIRO.AU Bill.Venables at CMIS.CSIRO.AU
Mon Oct 1 10:22:25 CEST 2001



>  -----Original Message-----
> From: 	Redding, Matthew [mailto:ReddinM at prose.dpi.qld.gov.au] 
> Sent:	Monday, 1 October 2001 3:45 PM
> To:	'r-help at lists.r-project.org'
> Subject:	[R] Grouped Regression
> 
> HI All R Gurus, 
> 
> Please reply directly to this email if you have a suggestion, as I am not
> currently an R-help member.
> 
> I would like to complete a non-linear regression with groups.  Having gone
> down the path of completing such
> things under GENSTAT, I would prefer to do it with R.
> 
> My example data set is...
> 
>  TimeTreat    ExtrVol2     SorbedT
>          ET1        39.5       384.5
>          ET1        78.7       305.2
>          ET1       118.0       261.5
>          ET1       157.6       233.7
>          ET1       197.5       206.3
>          ET2        37.5       422.1
>          ET2        74.8       405.5
>          ET2       112.0       382.6
>          ET2       149.7       361.8
>          ET2       187.4       347.7
> 
> To do a straightforward nls on the data, some code like this would be
> necessary...
> 
>
char<-read.table("O:/R_STH/L_TBA/P_IndSer/IntLive/Ilem_Research/Matt/Incubat
> ion/Desorption Isotherms/RWork/DPIAds.txt", header=TRUE)

Uh?  I suppose so.  Wouldn't it be simpler to put the data set in your
working directory?  Better still, do this step once and save the image from
that point on.

> library(NLS)
> test1<-nls(SorbedT~A + B/(1 + D*ExtrVol2), data=char,
start=list(A=200,B=80,
> D=0.003),trace=TRUE)
> summary(test1)
> 
> 
> Given this code, How do I alter it to group regression based on the factor
> TimeTreat, and provide a comparison of the 
> results?

Not enough people realise that the coefficients in a non-linear regression
may be parametrized as indexed quantities.

Here is one way of doing it.

> tst <- nls(SorbedT ~ A[TimeTreat] + B[TimeTreat]/(1 +
D[TimeTreat]*ExtrVol2),
+ data = dp1, start = list(A = c(200, 200), B = c(80,80), D = c(0.001,
0.001)),
+ trace = T)
82946.47 :  2e+02 2e+02 8e+01 8e+01 1e-03 1e-03 
82367.47 :  2.738746e+02 1.916630e+02 9.237563e+00 9.091811e+01 2.423017e-03
9.647061e-04 
79252.85 :  2.883669e+02 1.750317e+02 1.333030e+00 1.126316e+02 1.604157e-02
9.099994e-04 
78725.84 :   2.755148e+02  1.422234e+02  2.988133e+01  1.552895e+02
-5.021803e-02  .414247e-04 
77101.27 :   2.601821e+02  7.970087e+01 -8.884658e+01  2.362855e+02
2.189281e-02  7.813934e-04 
52402.35 :   2.172250e+02 -3.010763e+01  4.124670e+01  3.784320e+02
3.783349e-02  7.532221e-04 
46948.38 :   2.135563e+02 -4.034542e+01  4.767352e+01  3.917020e+02
1.571341e-02  7.530422e-04 
40955.72 :   2.053626e+02 -6.017882e+01  7.337312e+01  4.174103e+02
1.410275e-02  7.527228e-04 
31148.60 :   1.898579e+02 -9.736439e+01  1.218645e+02  4.656113e+02
1.342821e-02  7.522047e-04 
17465.68 :   1.626501e+02 -1.624220e+02  2.068017e+02  5.499460e+02
1.326065e-02  7.515317e-04 
4388.065 :   1.218128e+02 -2.599712e+02  3.342096e+02  6.764111e+02
1.324469e-02  7.509433e-04 
37.07059 :   8.097130e+01 -3.574638e+02  4.616163e+02  8.028200e+02
1.324616e-02  7.507585e-04 
37.07051 :   8.097170e+01 -3.574484e+02  4.616165e+02  8.028046e+02
1.324581e-02  7.508045e-04 
> coef(tst)
           A1            A2            B1            B2            D1 
 8.097170e+01 -3.574484e+02  4.616165e+02  8.028046e+02  1.324581e-02 
           D2 
 7.508045e-04 

A2 looks a bit sick if you were expecting a positive asymptote.

Bill Venables, 
CMIS, CSIRO Marine Laboratories, 
PO Box 120, Cleveland, Qld. 4163
Phone: +61 7 3826 7251  Fax: +61 7 3826 7304
<mailto: Bill.Venables at cmis.csiro.au>
http://www.cmis.csiro.au/bill.venables/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list