3
$\begingroup$

I am doing factorial experiments in R. I noticed when I use my variables as they are vs. coding them into -1,1, they are all very different. Here is my sample code. This dataset is adapted from the Experimentation for Improvement online course in Coursera.

#Example 1: Non-coded, as they are
runner1 <- data.frame(
  A = c(5, 10, 5, 10, 5, 10, 5, 10), 
  B = c(4, 4, 9, 9, 4, 4, 9, 9), 
  C = c(400, 400, 400, 400, 500, 500, 500, 500), 
  Y= c(39, 35, 40, 41, 40, 38, 41, 42))
runner1
runner1lm <- lm(runner1$Y ~ runner1$A*runner1$B*runner1$C)
runner1lm

#Example 2:Coded (-1, 1)

runner2 <- data.frame(
  A = c(-1, 1, -1, 1, -1, 1, -1, 1), 
  B = c(-1, -1, 1, 1, -1, -1, 1, 1),
  C = c(-1, -1, -1, -1, 1, 1, 1, 1), 
  Y= c(39, 35, 40, 41, 40, 38, 41, 42))
runner2
runner2lm <- lm(runner2$Y ~ runner2$A*runner2$B*runner2$C)
runner2lm

Which should I use? The coded or non-coded one? Also, how do I transform my entire dataset from coded to non-coded, if ever? Do I need a package for that? I was thinking about inputting coding parameters (like all 5 in column A will be -1, all 10 will be 1, so on and so forth like that) then letting R do its job and, boom!, I have a transformed data frame full of coded variables (except the dependent variable Y ofc). Basically I am asking if there's any code or package that can magically turn runner1 to runner2?

$\endgroup$
1
  • $\begingroup$ Look at the factor function and its relatives. $\endgroup$
    – whuber
    Jan 27 at 17:01

2 Answers 2

8
$\begingroup$

In the first model, you have a difference of 5 units between the two types of cases for predictors A and B, and 100 units for predictor C. In the second model, you've just scaled all the predictors to have a 2-unit difference between the types of cases.

The models aren't fundamentally different. They just use different scales for the predictor variables. It's like using distance as a predictor, but with one model measuring distance in millimeters and another in kilometers. The relative scales and coefficients are just inversely related.

If your predictor variables are fundamentally continuous, it's generally poor practice to convert them to categories as you seem to be doing in your question. That also removes the natural interpretation of regression coefficients as a change in outcome per unit change of a predictor value in its original scale. From that perspective, I see no need to transform data in a way needed to produce a model like runner2lm.

$\endgroup$
4
  • $\begingroup$ But what if I have mixed? Let's say I am working with two treatments. One is continuous the other is binary categorical (like with x vs. no x). Do I need to code them? $\endgroup$
    – Deus Sema
    Jan 28 at 2:40
  • $\begingroup$ @DeusSema let the software code the categorical predictor correctly (as whuber said, “look at the factor function and its relatives”) and keep the continuous predictor in its conventional, original units. In my experience that is the convention for publication. $\endgroup$
    – EdM
    Jan 28 at 3:54
  • $\begingroup$ @whuber I turned them into as.factor and they still give the different coefficient values as runner2lm $\endgroup$
    – Deus Sema
    Jan 28 at 9:44
  • 1
    $\begingroup$ @DeusSema factors are coded internally as 0/1 instead of -1/1 as in your question. Use the model.matrix() function to see the internal numeric coding. Furthermore, coefficients for predictors involved in interactions are for the situation when interacting predictors are at 0. So the displayed coefficient values will differ. Predictions from all models will still be identical. See this answer on centering continuous predictors (as R assumes in runner1lm and runner2lm). Comparisons of coding are easier to see if you use "+" instead of "*". $\endgroup$
    – EdM
    Jan 28 at 14:57
7
$\begingroup$

The models are identical.

One way to see this (although it's not a proof) is to compare their predictions and note that they are equal.

rbind(`First model`=predict(runner1lm), `Second model`=predict(runner2lm))
              1  2  3  4  5  6  7  8
First model  39 35 40 41 40 38 41 42
Second model 39 35 40 41 40 38 41 42

Of course the coefficients differ because you have used different methods to express the data with numbers.

We would expect the outputs to be interconvertible, though. You can figure out how by writing the model predictions as formulas. Let $A,B,C$ be the $\pm 1$ codes used in the second model. In this model, expanding the interactions gives

$$y = \beta_0 + \beta_A A + \beta_B B + \beta_C C + \beta_{AB} A\times B + \beta_{AC}A\times C + \beta_{BC} B\times C + \beta_{ABC}A\times B \times C.$$

There are eight coefficients to estimate. Their estimates, in the order given, are

$$(\beta_0, \beta_A, \ldots, \beta_{ABC}) = \left(39\frac{1}{2}, -\frac{1}{2}, 1\frac{1}{2}, \frac{3}{4}, 1, \frac{1}{4}, -\frac{1}{4}, -\frac{1}{4}\right),$$

as you saw when you printed runner2lm at the end of your code.

For instance, in the first record where $(A,B,C)=(-1,-1,-1),$ the prediction is

$$y = 39\frac{1}{2} -\frac{1}{2}(-1) + \cdots + -\frac{1}{4}(-1)(-1)(-1) = 39.$$

The first model changes how $A,$ $B,$ and $C$ are expressed. Let's call these variables $a,b,$ and $c.$ The relationships are

$$\left\{\begin{aligned} a &= \frac{15}{2} + \frac{5}{2}A;&\quad A &= (2a - 15)/5\\ b &= \frac{13}{2} + \frac{5}{2}B;&\quad B &= (2b - 13)/5\\ c &= 450 + 50C &\quad C &= (c - 450)/50. \end{aligned}\right.$$

Plugging these expressions for $A,B,C$ into the model gives

$$y = \beta_0 + \beta_A(2a-15)/5 + \cdots + \beta_{ABC}\left((2a-15)/5\right)\left((2b-13)/5\right)\left((c-450)/50\right).$$

After expanding all these products and doing quite a bit of algebra, you find the corresponding coefficients for the first model. For instance, the coefficient of $abc$ is $$\beta_{ABC}(2/5)(2/5)(1/50) = -4/(4\times 5\times 5\times 50) = 1/1250 = 0.0008.$$

You can inspect this directly:

tail(coefficients(runner1lm), 1)
runner1$A:runner1$B:runner1$C 
                       -8e-04

-8e-04 is computerese for $8\times 10^{-4} = 0.0008,$ agreeing with the calculation.

It is a worthwhile exercise to check some of the other coefficients in the same way: inspect the coefficient in one model; do the algebra to work out what it must be in the other model; and then confirm that through inspection. Repeat this exercise until you feel you really understand the relationships between these two forms of the same model.

If you get stuck, simplify your example. Start with a model of the form $y = \beta_0 + \beta_A A,$ with one explanatory variable. Then move to a model with two variables and an interaction term. After that you should have the hang of it.

Moral

When you report model estimates, make sure to explain how you chose numbers to express all the variables involved.

Some methods of numerical expression ("coding") lend themselves better to interpretation than others, but any method that creates a mathematically predictable one-to-one correspondence, as examined above, will work.

Computing bonus

Because the relationships between these forms of numerical expression can be computed algebraically, they permit easy, efficient transformation between the data frames. For instance, the expression

15/2 + (5/2)*runner2$A

converts the A variable in runner2 (which I named $A$ in the formulas above) to the A variable in runner1 (which I named $a$). Thus, you don't need to create any translation tables or use any special packages: addition and multiplication take care of everything.

$\endgroup$
3
  • $\begingroup$ What is the convention in publications? Coded or non-coded? $\endgroup$
    – Deus Sema
    Jan 28 at 2:37
  • $\begingroup$ The convention is what I wrote under "Moral." $\endgroup$
    – whuber
    Jan 28 at 14:04
  • $\begingroup$ Sorry the moral part sounds very vague for me. I was expecting a direct answer whether it should be coded or non-coded as per journal and publication convention. $\endgroup$
    – Deus Sema
    Feb 1 at 14:41

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.