1
$\begingroup$

Data

I'm constructing linear models using the fish dataset from this data with the following variables:

  • res_var: response variable Weight
  • exp_var1: explanatory variable Length
  • exp_var2: explanatory variable Species
library(dplyr)

fish <- read_fst("Fish.csv")

fish <- fish %>%
mutate(Length = (Length1 + Length2 + Length3) / 3)

Model 1

I am trying to understand the difference between the coefficients of the following two models. The first is a model with two explanatory variables and their interactions:

lm(Weight ~ Length + Species + Species:Length + 0, data = fish)

Output:

Coefficients:
             Length         SpeciesBream         SpeciesPerch          SpeciesPike         SpeciesRoach  Length:SpeciesPerch   Length:SpeciesPike  Length:SpeciesRoach  
             50.869            -1107.791             -640.845            -1560.734             -342.989              -13.979               -0.836              -28.994 

Model 2

The second model includes just one explanatory variable and its interaction with another explanatory variable:

lm(Weight ~ Species + Species:Length + 0, data = fish)

Output:

Coefficients:
       SpeciesBream         SpeciesPerch          SpeciesPike         SpeciesRoach  SpeciesBream:Length  SpeciesPerch:Length   SpeciesPike:Length  SpeciesRoach:Length  
           -1107.79              -640.85             -1560.73              -342.99                50.87                36.89                50.03                21.88  

To understand the second model above, I first filtered the fish dataset for each Species:

bream <- fish %>% 
  filter(Species == "Bream")

perch <- fish %>% 
  filter(Species == "Perch")

pike <- fish %>% 
  filter(Species == "Pike")

roach <- fish %>% 
  filter(Species == "Roach")

Then, I constructed a simple linear model for each subset dataset above:

mdl_bream <- lm(Weight ~ Length, data = bream)

mdl_perch <- lm(Weight ~ Length, data = perch)

mdl_pike <- lm(Weight ~ Length, data = pike)

mdl_roach <- lm(Weight ~ Length, data = roach)

I then looked at their coefficients.

coefficients(mdl_bream)

coefficients(mdl_perch)

coefficients(mdl_pike)

coefficients(mdl_roach)

Output:

(Intercept)      Length 
-1107.79052    50.86892 
(Intercept)      Length 
 -640.84506    36.89006 
(Intercept)      Length 
-1560.73352    50.03289 
(Intercept)      Length 
 -342.98917    21.87535 

These correspond to the second model I mentioned: lm(Weight ~ Length + Species + Species:Length + 0, data = fish).

Can anyone help me understand how the first model's coefficients should be interpreted?

$\endgroup$
4
  • $\begingroup$ I added the interaction tag. You can browse it for answers to your question (I don't hve much time right now). I also removed the R tag as this doesn't really have anything to do with R. $\endgroup$
    – Peter Flom
    Apr 1 at 23:21
  • $\begingroup$ @PeterFlom Thanks for your help. I added the R tag, thinking syntax might be important here. $\endgroup$ Apr 1 at 23:27
  • 1
    $\begingroup$ @PeterFlom I have added it back because I think it is difficult to understand what they are fitting here without knowing R. Otherwise I would agree that the question is largely about fitting models over programming. $\endgroup$ Apr 1 at 23:34
  • $\begingroup$ An aside - in R * and : produce different results, which can add to the confusion of building models 😜 $\endgroup$
    – N Brouwer
    Apr 3 at 4:34

2 Answers 2

1
$\begingroup$

You should almost always include the main effects when fitting an interaction. Some call this the "hierarchical principle". Logically, it is odd to only consider an interaction independent of the main effects, as if the main effects don't exist but only the interaction between the two variables does. To me a categorical-by-continuous interaction is more sensible anyway. It represents the adjusted slope based off the grouping variable. And surely here this would make sense given we are talking about fish, which undoubtedly have differences in weight and length by species.

Edit 1

Here is an R example with the likely familiar iris dataset, which shows how categorical-by-numeric variable interactions work. Here I first fit the data and inspect the coefficients:

#### Fit Model ####
fit <- lm(
  Petal.Length ~ Species * Petal.Width,
  iris
)

coefs <- coef(fit)
coefs

The coefficients seem to indicate that the setosa flowers are the reference group and the others are compared against it, with the slope and intercept increasing based off their coefficients. For example, virginica seems to have a larger mean (by it's larger coefficient by itself) and a slightly more positive slope (indicated by the interaction term):

                  (Intercept)             Speciesversicolor              Speciesvirginica                   Petal.Width 
                    1.3275634                     0.4537120                     2.9130892                     0.5464903 
Speciesversicolor:Petal.Width  Speciesvirginica:Petal.Width 
                    1.3228344                     0.1007691 

This becomes more clear when you simply plot the model terms directly. Here I plot the regression fit we just used, but explicitly code the regression lines to fit the intercepts and slopes by each group. I have annotated the geom_abline code to show what is being done here:

#### Plot Model ####
iris %>% 
  ggplot(aes(x=Petal.Width,
             y=Petal.Length,
             color=Species))+
  geom_point(size = 5,
             alpha = .4)+
  geom_abline(
    intercept = coefs["(Intercept)"], # conditional mean of baseline group "setosa"
    slope = coefs["Petal.Width"], # and it's specific slope
    color = "red",
    linewidth = 1
  )+
  geom_abline(
    intercept = coefs["(Intercept)"] + coefs["Speciesversicolor"], # conditional mean adjustement
    slope = coefs["Petal.Width"] + coefs["Speciesversicolor:Petal.Width"], # slope adjustment
    color = "green",
    linewidth = 1
  )+
  geom_abline(
    intercept = coefs["(Intercept)"] + coefs["Speciesvirginica"], # conditional mean adjustment
    slope = coefs["Petal.Width"] + coefs["Speciesvirginica:Petal.Width"], # slope adjustment
    color = "blue",
    linewidth = 1
  )+
  labs(x="Petal Width",
       y="Petal Length",
       title="Petal Dimensions by Species")+
  theme_classic()

And you can see now the regression lines change based off each group:

enter image description here

Edit 2

In the comments you also noted whether or not one should fit models which use subsets over the full dataset. The first problem with this is that you increase the likelihood of a Type I error by fitting multiple models. Second, I would think a model which estimates direct comparisons of groups would be more meaningful than a handful of separate ones. Fitting them separately is weird. It's like saying "the model for fish measures in sea bass is a completely different thing from the model for fish measures in trout." This is obviously not true and wouldn't really be meaningful in my opinion. Some similar questions on this topic can be found here and here, where the answers largely echo my own sentiments.

Edit 3

GG's answer is correct regarding the specification of the model syntax (as I note in the comments of his answer, I missed this on my first read). However, the rest of what I say here is still true (particularly the part about subsetting your data to fit separate models).

$\endgroup$
5
  • $\begingroup$ Thanks for sharing the hyperlink to the "hierarchical principle." I apologize if that's a naive question, but how is the slope adjusted based on a categorical variable? Edit: are you suggesting Model 1 is more appropriate? $\endgroup$ Apr 2 at 0:06
  • $\begingroup$ See my edit, where I use an R example using the iris dataset. $\endgroup$ Apr 2 at 0:19
  • $\begingroup$ Thanks so much! What's the advantage of using fit over the following setosa <- iris %>% filter(Species == "setosa") lm(Petal.Length ~ Petal.Width, data = setosa), etc. $\endgroup$ Apr 2 at 0:35
  • $\begingroup$ I have edited my answer again to clarify that point. $\endgroup$ Apr 2 at 0:51
  • 1
    $\begingroup$ These make a lot of sense. Thanks a lot for sharing the additional resources. From now on, I'll avoid splitting the data or excluding the main effects. Thanks again! $\endgroup$ Apr 2 at 0:51
0
$\begingroup$

Consider the 5 models below. fm1 is Model 1 from the question. fm2 is just a different way of writing it. fm3 uses / notation and is similar to fm4 but is a single model. fm4 is equivalent to running separate models for each level of Species and gives a list of models as the result. fm5 is model 2 from the question.

Now

  • fm1 and fm2 have the same coefs
  • fm3 and fm4 have the same coefs
  • all of them have the same fitted values, residuals and deviances (aka residual sum of squares) . Note that if the coefs are the same then the fitted values are necessarily the same. Also if the fitted values are the same then since the residuals equal the dependent variable minus the fitted values then the residuals are the same. Also if the residuals are the same then the residual sum of squares, being a function of the residuals, must be the same.

so all of them give the same predictions. Thus you can use any one of them depending on whichever you find easiest to think about.

fm1 <- lm(Weight ~ Length + Species + Species:Length + 0, fish)
fm2 <-  lm(Weight ~ Length * Species + 0, fish)
all.equal(coef(fm1), coef(fm2)) # TRUE

fm3 <- lm(Weight ~ Species / Length + 0, fish)
fm4 <- nlme::lmList(Weight ~ Length | Species, fish)
all.equal(unname(coef(fm3)), unname(unlist(coef(fm4)))) # TRUE

fm5 <- lm(Weight ~ Species + Species:Length + 0, fish)

all.equal(fitted(fm1), fitted(fm2)) # TRUE
all.equal(fitted(fm1), fitted(fm3)) # TRUE
all.equal(fitted(fm1), fitted(fm5)) # TRUE

A subtle point regarding R lm formulas as they relate to Model 2 of the question, fm5 here, is that it may appear that the main effect of Length was omitted but in actuality it was not. Length can be expressed as a linear combination of the columns of the model matrix of Species:Length + 0 as this shows:

deviance(lm(Length ~ Species:Length + 0, fish)) # 6.825225e-27

Note

To create fish

library(dplyr)

u <- "https://raw.githubusercontent.com/Ankit152/Fish-Market/main/Fish.csv"
fish <- u |>
  read.csv() |>
  mutate(Length = (Length1 + Length2 + Length3) / 3)
$\endgroup$
5
  • $\begingroup$ To whoever downvoted this the results can easily be reproduced by running the code in the Note and then the code in the body of the answer if you don't believe it. It shows everything claimed is true. $\endgroup$ Apr 3 at 0:24
  • $\begingroup$ You're right about the code (I somehow missed that they simply replaced * with the +/: combo), but I still think I have an issue with the following statement: Thus you can use any one of them depending on whichever you find easiest to think about.. You can't just use any of them willy nilly. For the reasons I explained in my own answer, this can be a problematic practice if we fit the data to subsets only (which isn't addressed in your answer), particularly if you are trying to match a theoretical model to a mathematical model. $\endgroup$ Apr 3 at 1:53
  • $\begingroup$ You can use any of them because they all give the same predictions and only differ in parameterization. Regarding the subsets fm4 produces all subset models in a list and fm3 gives the exact same parameters as fm3 but using a single lm model and we further show that fm3 and hence fm4 give the same predictions as all the others. $\endgroup$ Apr 3 at 2:08
  • $\begingroup$ The coefficients will be the same, but not everything else. You are forgetting about model differences in SE, RSE, $R^2$, and statistical power. If you care about the $p$ values, then this is also problematic to do. When using the conventional alpha cutoffs, this increases your probability of committing a Type I error from 5% to 20% (if using the four subsetted models). It also tests different hypotheses (such as the null of the model of a full fit vs a subsetted one). I also don't know why you would want a model that is much more limited in scope for future uses (if prediction is the goal). $\endgroup$ Apr 3 at 2:40
  • $\begingroup$ Some (not all) have different coefs so clearly we can't compare the std errs but the key point is that the models are all equivalent from the viewpoint of giving the same predictions. Also contrary to the claim, the R^2 and sigma are the same. R^2 equals cor(Weight, fitted(fm1))^2 and since Weight and fitted(fm1) are the same for fm1, fm2, fm3 and fm5 they all have the same R^2. summary(fm1)$r.squared equals 0.9798636 and the same for fm2, fm3 and fm5. sigma(fm1) equals 79.46825 and the same for fm2, fm3 and fm5. fm4 too if you use the pool=TRUE argument to lmList. $\endgroup$ Apr 3 at 7:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.