From the graphs it appears that the relationship between Body weight and Heart weight is linear, so we can try to fit a linear model. One way is to use the method of least squares regression:
In R we use the command lm(y~x) to find the least squares regression line. To do it for the cats data and for the females only we could use the command lm(Hwt[Sex=="F"]~Bwt[Sex=="F"]) or alternatively lm(Hwt~Bwt, subset= Sex=="F").
If we run fit <-lm(y~x), then "fit" is a special, rather complicated data structure in R. Such objects often have methods already implemented for them. For example an object of class "lm" has the following methods:
• summary() prints useful information.
• coef() gives the coefficients of the model (b0, b1,..).
• resid() finds the residuals
• fitted() finds the fitted values
• deviance() finds the residual sum of squares
• anova() finds the sequential analysis of variance table
• predict() predicts the mean response for new observations
• plot() draws a series of diagnostic plots
Note that in this problem it might make good sense to fit a no-intercept mode of the form y=b1x (why?) This can be done using the command lm(Hwt~Bwt-1, subset= Sex=="F")
How about fitting a linear model for both sexes? One way to do this is to just use the command above with Sex=="M" instead of Sex=="F" to get the fit for the males. Another is to use the variable "Sex" as an additional variable, that is fitting the model Hwt~Bwt+Sex. It is easy to see, though, that this fits two parallel lines, which is clearly wrong here. Instead we need to include an interaction term. We can use the command lm(Hwt~Bwt+Sex+Sex*Bwt). The command lm(Hwt~Sex*Bwt) does the same thing.
What are these lines if we write them down for males and females separatey?
The same model can be fitted explicitly, and withoud the need for any arithmetic, using the command lm(Hwt~Sex/Bwt-1)
In cats.fun with which=4 we draw the least squares regression lines in the scatterplot.
How good a model is this? Specifically, is the relationship really linear? A first look at this is given by the residuals vs. fits plot. It is drawn as part of a series of plots by the plot() command, and it indicates a good model here.
Using the method of least squares to find the equation is always "legal" but not always a good idea. lm.exa(1) shows some cases where the least squares regression line is clearly not a good model for the data.
All the other parts of the output of summary() have the additional condition ei~N(0,s2). Notice that this has two parts:
1) The residuals have a normal distribution
2) The standard deviations are all the same (homoscatasticity)
We can check these as follows:
1) Draw the normal probability plot of the residuals (It is already part of the plot() routine)
2) Check the residuals vs. fits plot again and look for a change in variance from left to right. lm.exa(2) shows an example with non-constant variance.
For the cats data the normal assumption appears to be fine but there is a slight hint of unequal variance.
One way to try and deal with the problem of unequal variance is to use a transformation. In our case we really are interested in testing Hwt~Bwt. If this is so then Hwt/Bwt should be a constant. So we will now consider the model log(Hwt/Bwt)=b0+b1log(Bwt)+e
The scatterplot of log(Hwt/Bwt) vs. log(Bwt) shows no indication of a relationship.
fit <- lm(log(Hwt/Bwt)~Sex/log(Bwt),data=cats)
summary(fit) has an F-statistic of 1.46 on 3 and 140 df with a p-value of 0.229, which indicates that the total contribution of the model beyond the constant is negligible, suggesting that in both sexes heart weight is proportional to body weight and that the constant is proportionality is the same for both sexes.
Using the F-statistic for this question is not the same as looking at the individual t-tests. Why not?
A check of the diagnostic plots shows that the log transformation was succesful in removing the heteroscatasticity.