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)) obs exp.x -0.04791535 -0.04791535
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)) 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!