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]-mean(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
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]-mean(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 >
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 >