2
$\begingroup$

I have data that consists of a grouping variable $\tt{grp}$, a predictor $\tt{x}$, and a response $\tt{y}$. There are three groups. I want to fit to this data the model $y = \alpha_i + (\beta + \gamma_i)x + \epsilon_i$, where $i = 1, 2, 3$ is the group number, $\sum_{i = 1}^3 \gamma_i = 0$, and $\epsilon_i$ is an error term. Each group has its own intercept and slope; the deviations $\gamma_1, \gamma_2, \gamma_3$ of the group slopes from some baseline $\beta$ sum to zero. Is there a way to fit this model using R's lm function? Below is some code that generates data like mine.

replicate(
  3,
  {
    tibble(
      x = seq(1, 2, by = 0.1),
      y = rnorm(1, sd = 0.1) + rnorm(1, sd = 1) * x + rnorm(11, sd = 0.05)
    )
  },
  simplify = FALSE
) %>%
  bind_rows(.id = "grp") %>%
  mutate(grp = as_factor(grp)) ->
  tbl
$\endgroup$
2
$\begingroup$

If I understand you correctly, you want to apply "sum to zero contrasts". That is not constrained regression but just a reparameterization of the model fit with the default treatment contrasts. The model fit (RSS/likelihood, predictions, etc.) is exactly the same.

Note that only two of the $\gamma_i$ will be estimated by the model since the third one is implicit from the sum to zero. Also, the following applies these to the intercept, too. But you can easily calculate the $\alpha_i$ (see below).

fit <- lm(y ~ x * grp, data = tbl, contrasts = list(grp = contr.sum))    
coef(fit)
#(Intercept)           x        grp1        grp2      x:grp1      x:grp2 
#-0.03465095  0.56939954 -0.12084281  0.08225554  0.65901278 -0.72977953

Illustrate how to interpret the output:

predict(fit, newdata = data.frame(x = 0, grp = factor(1:3)))
#           1            2            3 
#-0.155493755  0.047604589  0.003936318
-0.03465 + -0.12084
#[1] -0.15549
-0.03465 + 0.08226
#[1] 0.04761
-0.03465 + (0 - 0.08226 - -0.12084)
#[1] 0.00393

predict(fit, newdata = data.frame(x = 1, grp = factor(1:3)))
#         1          2          3 
# 1.0729186 -0.1127754  0.6441026
-0.03465 + -0.12084 + 0.56940 + 0.65901
#[1] 1.07292
-0.03465 + 0.08226 + 0.56940 + -0.72978
#[1] -0.11277
-0.03465 + (0 - 0.08226 - -0.12084) + 0.56940 + (0 - 0.65901 - -0.72978)
#[1] 0.6441 
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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