# Calculating regression coefficients from the covariance matrix involving categorical variables

## 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!

• Is your categorical variable ordinal one and you want to think there is a normally distributed underlying variable "behind" it? - so that you are speaking of polyserial correlation? Dec 27, 2018 at 12:45
• Yes, that's right. Dec 27, 2018 at 12:53
• The relevant covariance matrix is that of the model matrix. There's no reason to suppose this would be identical to that produced by the polyserial correlation. You might want to inspect the model matrix and its covariance to identify the differences.
– whuber
Dec 27, 2018 at 16:09
• @whuber, thank you for the insight. Using the model matrix does yield the same results. I'm still studying where the differences are, why they happen, and how to work this out in my actual case. I'll add my findings ASAP. Dec 28, 2018 at 10:36

## 1 Answer

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