Calculating regression coefficients from the covariance matrix involving categorical variables

Question

Prelude

Let us start by generating some data:

set.seed(271317)
x <- rnorm(n = 1e3)  # independent variable
y <- rnorm(n = 1e3)  # response variable (indep. from X; shouldn't matter here)

I am interested in estimating $\beta$ in the following model:

$$ \hat y = \alpha + \beta x $$

One way to do is is using the covariance matrix between $Y$ and $Y$. In this simple case:

beta_x <- cov(y, x) / var(x)

In R, we can also use lm() to fit a linear model and retrieve the coefficients.

reg_yx <- lm(y ~ x)

And the results are pretty similar, as expected:

> print(c(obs = beta_x, exp = reg_yx$coef[2]))
        obs       exp.x 
-0.04791535 -0.04791535 

---

Problem

Let us now introduce a categorical ordinal variable $W$:

prob_w <- .2
w <- as.factor(rbinom(n = 1e3, size = 1, prob = prob_w))

I understand that calculating the slope for this regression using the covariance matrix is still possible, as long as one uses the polyserial correlation (an inferred correlation) instead of good ol' Pearson correlation. From the polyserial correlation I get the covariance:

var_w  <- prob_w * (1 - prob_w)
cor_yw <- polycor::polyserial(y, w)  # approx. 0.01969907
cov_yw <- cor_yw * sqrt(var_w) * sd(y)

Then $\beta_W$ is calculated as:

beta_w <- cov_yw / var_w

As a benchmark, I am using a linear regression through lm() again, which automatically dummy codes $W$:

reg_yw <- lm(y ~ w)

However, now the betas don't match anymore. In some cases, they differ by a whole order of magnitude, so I don't think sample variability is to blame here.

print(c(obs = beta_w, exp = reg_yw$coef[2]))
       obs     exp.w1 
0.04804846 0.03344181 

I bet I am missing something simple, but I've been stuck on this problem for way too long to be able to see it anymore. Help!

Answer 1

Thanks to the comment from @whuber, the solution requires using the model matrix instead of the original vector w of factors.

set.seed(271317)
x <- rnorm(n = 1e3)  # generated so that y and w also match the OP
y <- rnorm(n = 1e3)
prob_w <- .2
w <- as.factor(rbinom(n = 1e3, size = 1, prob = prob_w))

reg_yw <- lm(y ~ w)

model_w <- model.matrix(y ~ w)
cov_yw_model <- cov(model_w, y)[-1]  # removes the intercept
var_w_model <- var(model_w)[-1, -1]  # removes the intercept
beta_w_model <- cov_yw_model / var(model_w)[-1, -1]

The following shows the results are identical.

print(rbind(obs = beta_w_model, exp = reg_yw$coef[-1]))
            w1
obs 0.03344181
exp 0.03344181

---

For what it's worth, the solution above also works for a regression on a combination of continuous, binomial and polinomial variables. We just have to remember to use solve() to calculate the slopes. Here's a quick example:

set.seed(271317)
n <- 100

y <- rnorm(n)
x <- rnorm(n)
w <- as.factor(rbinom(n, size = 1, prob = .5))
z <- as.factor(rbinom(n, size = 2, prob = .5))

reg <- lm(y ~ x + w + z)

model <- model.matrix(y ~ x + w + z)
cov_model   <- var(model)[-1, -1]
cov_model_y <- cov(model, y)[-1]
betas <- solve(cov_model, cov_model_y)
print(rbind(obs = betas, exp = reg$coef[-1]))

The final line yields

             x         w1        z1         z2
obs -0.0249083 0.03175005 0.0634217 -0.1570229
exp -0.0249083 0.03175005 0.0634217 -0.1570229