[R] Arrays of functions or results of functions.

[Berton Gunter]

If I understand you correctly, I would say that this is a standard analysis
of covariance problem that you are approaching incorrectly. You should not
be testing for "equal variances," IMO. Instead,

1. Combine all your data into 3 columnsS x, y, and group= subset.

2. Model.1: y ~ x
3. Model.2: y ~ x*group.  This gives a separate slope and intercept for each
group.
4. anova(Model.1, Model.2)  to compare.

A better approach would be to use lmList() from the nlme package. Better yet
would be model the data as a mixed effect model via lme(), which assumes
that slopes and intercepts vary randomly among your subsets about their
corresponfding central values. This gets you shrinkage, which is highly
desirable. See Bates and Pinheiro "MIXED EFFECTS MODELS IN S AND S-PLUS" for
details.

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
 
 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Gabor 
> Grothendieck
> Sent: Monday, March 27, 2006 9:48 AM
> To: Gary Mallard
> Cc: r-help at stat.math.ethz.ch
> Subject: Re: [R] Arrays of functions or results of functions.
> 
> Is your question how to run a regression with separate slopes
> and then with one slope and then to complare them?  If that
> is it here is an example:
> 
> # test data - ind is factor which defines the groups
> set.seed(1)
> y1 <- 10 + 20 * seq(100) + rnorm(100)
> y2 <- -200 + 35 * seq(100) + rnorm(100)
> yy <- pmax(y1, y2)
> ind <- factor(y1 > y2)
> 
> # lm.N has N slopes
> lm.1 <- lm(yy ~ ind + seq(100) - 1)
> lm.2 <- lm(yy ~ ind/seq(100)-1)
> lm.1
> lm.2
> anova(lm.1, lm.2)
> 
> On 3/27/06, Gary Mallard <gary.mallard at nist.gov> wrote:
> > The general problem I am trying to solve is to determine if 
> a series of
> > subsets of data can be described with a single regression 
> slope.  This
> > involves fitting the data to each subset, calculating a joint slope
> > followed by F tests to show that the variances are equal 
> the final slope is
> > valid.
> > The data for is characterized by a parameter PC (for the 4 
> - in this case)
> > sets and a dependent variable RI and an independent variable ROH.
> > The data are contained in a variable "joint".
> > The joint has been attached and has RI, ROH and PC for each element.
> > The following gives the initial results:
> > Mline<-lm(RI[PC==1]~ROH[PC==1])
> > Eline<-lm(RI[PC==2]~ROH[PC==2])
> > Iline<-lm(RI[PC==3]~ROH[PC==3])
> > Pline<-lm(RI[PC==4]~ROH[PC==4])
> >
> >
> > joint_reduced <- joint;
> > for(i in 1:4) {
> >        
> joint_reduced$RI[joint_reduced$PC==i]<-joint$RI[joint$PC==i]-m
> ean(joint$RI[joint$PC==1]}
> > AllLine<-lm(joint_reduced$RI~joint_reduced$ROH);
> >
> > Now the statistics from AllLine can be compared with each 
> of the individual
> > statistics.
> >
> > NOW THE QUESTION:
> >  From a lot of point of view it would be useful to have a parameter
> > generated by
> > for (i in 1:4){         Xline[i]=lm(RI[PC==i]~ROH[PC==i])}
> > And now all of the work of comparison can be done with 
> calls to Xline[i]
> > rather than having to  work individually with Mline, Eline, etc.
> >
> > This appears to be impossible.  The constructor for Xline[i] is not
> > automatic (as it is for Mline, etc) noted above.  I cannot 
> determine how to
> > construct the Xline[i] object so that this kind of process can be
> > generalized.  Is it possible?  Is there another way to set 
> us such tests of
> > multiple line linearity that is already in a package?
> >
> > Comments or pointers would be appreciated.
> > Gary
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> >
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>