Hi everyone :] first time posting here. I have a question about linear models. Suppose I have two models of the form:
$y_1 = \beta_1x + error_1$
$y_2 = \beta_2x + \beta_3y_1 + error_2$
So graphically, the system looks like:
And I’m interested in estimating the coefficient $\beta_3$, which is fixed. $\beta_1$ and $\beta_2$, meanwhile, are random variables.
The natural thing to do would be to fit the second model. But I don’t have access to any of the ’x’s or ’y’s -- instead, I have lots and lots of estimates of $\beta_1$s from the first model, as well as $\beta_4$s from the model:
- $y_2 = \beta_4x + error_3$
When I regress $\beta_4$ ~ $\beta_1$, ie fit:
- $\hat{\beta_4} = \beta_5\hat{\beta_1} + error_4$
the estimate I get for the linear coefficient describing their relation, $\hat{\beta_5}$, appears to converge upon $\beta_3$ at high n. This seems to hold in simulation across lots of different values for different coefficients. In R:
library(MASS)
b3 = round(rnorm(1,0,2),2)
# b2 = round(rnorm(1,0,2),2)
b2 = NA
sim <- function(b2 = NA, b3 = NA){
n = 1E3
x <- rnorm(n, sd = rexp(1, 1))
b1 <- rnorm(1, sd = rexp(1, 1))
if(is.na(b2)){b2 <- rnorm(1, sd = rexp(1, 1))}
if(is.na(b3)){b3 <- rnorm(1, sd = rexp(1, 1))}
y1 <- x*b1 + rnorm(n, sd = rexp(1, 1))
y2 <- x*b2 + b3*y1 + rnorm(n, sd = rexp(1, 1))
fit1 <- lm(y1 ~ x)
fit2 <- lm(y2 ~ x)
c(b1 = fit1$coefficients["x"], b4 = fit2$coefficients["x"])
}
fits <- t(replicate(1E4, sim(b2, b3)))
plot(fits)
coef_fit <- rlm(fits[,"b4.x"] ~ fits[,"b1.x"])
abline(coef_fit, col = 2, lty = 2)
text(paste0(paste0(rep(" ", 100), collapse = ""), "true value of b3 = ", b3, ", estimated slope = ", round(coef_fit$coefficients["fits[, \"b1.x\"]"], 3)),
x = par("usr")[1], y = par("usr")[3], pos = 3)
text(paste0(paste0(rep(" ", 100), collapse = ""), "true value of b2 = ", ifelse(is.na(b2), "Exp. 0", b2), ", estimated intercept = ", round(coef_fit$coefficients["(Intercept)"], 3)),
x = par("usr")[1], y = par("usr")[4], pos = 1)
(I’m using rlm() , whatever it’s doing — presumably using a Student’s t in place of the normal in the error term, which seems poetic, if nothing else, since the coefficients themselves ought be t-distributed — and also more practically, cos the coefficient estimates are noisy, depending as they do on variables drawn from a variance mixture of normals)
The broader intuition here, was that some of the effect of $\beta_3$ will manifest as extra effect captured by $\hat{\beta_4}$ in proportion to $\beta_1$, and so examining the distribution of $\hat{\beta_4}$s to $\hat{\beta_1}$ will tell you something about $\beta_3$ (and $\beta_2$). Specifically, if we substitute 1) into 2) we get:
$y_2 = \beta_2x + \beta_3(\beta_1x + error_1) + error_2$
$y_2 = \beta_2x + \beta_3\beta_1x + \beta_3error_1 + error_2$
$y_2 = (\beta_2+ \beta_3\beta_1)x + error_{heteroskedastic}$
and estimate equation 3), which assumes constant residual variance, to get an estimate of $\beta_4$, which multiplies x to get $y_2$. What else multiplies x to get $y_2$? Why, $\beta_2 + \beta_3\beta_1$ in equation 7). So if we set $E(\beta_4) = \beta_2 + \beta_3\beta_1$, regressing $\beta_4 \sim \beta_1$, the slope yields an estimate of $\beta_3$ and the intercept an estimate of $\beta_2$. As such, a lingering issue might entail the failure of 3) to accommodate the increasing residual variance with x, though apparently that doesn't bias our estimates of $\hat{\beta_4}$.
But I’m having trouble demonstrating more rigorously why this is the case -- my initial foray was to try to describe the joint distribution of all variables, and then work out the conditional distribution of $\beta_3$ given $\beta_1$ and $\beta_4$, themselves given x, $y_1$, and $y_2$, which are hidden from me. But in this scenario, $\beta_3$ is fixed, though I guess I could try treating it as a random variable in itself. Also, things like $x\beta_1$ are a product of normals here, which aren’t very nice to work with. And I tried working through things at the level of implied covariance patterns (e.g. cov(X, $Y_1$) = $\beta_1$var(x)... a bit of algebra later and $\beta_3 = \frac{cov(x,y_2)}{\beta_1var(x)} - \frac{\beta_2}{\beta_1}$), as well as formulae for partial correlations etc. but couldn’t get anything to stick, and also I don’t actually know $\beta_2$ and $\beta_1$, just $\hat{\beta_1}$).
I tried to look at the coverage properties of this ‘estimator’ and they’re not quite there, but not terrible either:
sim2 <- function(n){
b3 = round(rnorm(1,0,2),2)
fits <- t(replicate(n, sim(b3 = b3)))
coef_fit <- rlm(fits[,"b4.x"] ~ 0 + fits[,"b1.x"], maxit = 100)
in95CI <- abs(summary(coef_fit)$coefficients[1] - b3) < qt(0.975, df = n-1)*summary(coef_fit)$coefficients[2]
tscore <- (summary(coef_fit)$coefficients[1] - b3) / summary(coef_fit)$coefficients[2]
return(c(in95CI = in95CI, tscore = tscore))
}
n = 1E2
coverage_sims <- as.data.frame(t(replicate(500, sim2(n))))
mean(coverage_sims$in95CI)
hist(coverage_sims$tscore, breaks = seq(min(coverage_sims$tscore), max(coverage_sims$tscore), length.out = 50))
hist(pt(coverage_sims$tscore, df = n-1), breaks = seq(0,1,length.out = 20))
(for a pointier take, the 95% CI looks to cover the true value of $\beta_3$ about 79% of the time)
In any case -- is this a valid inferential procedure? If it is, can anyone point me to or provide a more formal derivation? If it’s not, can someone under what circumstances it would behave or misbehave?
(I'd also note that in my true use-case the 'x's would not be iid, but I'd have an rough estimate of their similarity / distance, so I'd probably need to use a GP for $\beta_4 \sim \beta_1$. Didn't want to complicate things too much in this question, though)
EDIT: also, I suspect similar reasoning can be used to extend this to a multiple regression framework, where:
$\vec{y} \sim MVN(\vec{\mu}, \Sigma)$
&
$z = \sum_{i=1}^{n}\beta_iy_i + error$
A quick numerical demonstration:
rlkj <- function (K, eta = 1) {
alpha <- eta + (K - 2)/2
r12 <- 2 * rbeta(1, alpha, alpha) - 1
R <- matrix(0, K, K)
R[1, 1] <- 1
R[1, 2] <- r12
R[2, 2] <- sqrt(1 - r12^2)
if (K > 2)
for (m in 2:(K - 1)) {
alpha <- alpha - 0.5
y <- rbeta(1, m/2, alpha)
z <- rnorm(m, 0, 1)
z <- z/sqrt(crossprod(z)[1])
R[1:m, m + 1] <- sqrt(y) * z
R[m + 1, m + 1] <- sqrt(1 - y)
}
return(crossprod(R))
}
n_y <- 20
bs_notfitted = round(rnorm(n = n_y,0,2),2)
corr_mat_ys <- rlkj(n_y, 1)
sim3 <- function(bs_notfitted, corr_mat_ys){
n_obs = 1E3
n_y <- length(bs_notfitted)
x <- rnorm(n_obs, sd = rexp(1, 1))
bs_fitted <- rnorm(n_y, sd = rexp(1, 1))
ys_exp <- crossprod(t(x), bs_fitted)
sd_ys <- rexp(n_y, 1)
ys <- t(sapply(1:n_obs, function(obs)
rmvnorm(n = 1, mean = ys_exp[obs,], sigma = diag(sd_ys) %*% corr_mat_ys %*% diag(sd_ys))
))
z <- x*b1 + ys %*% bs_notfitted + rnorm(n_obs, sd = rexp(1, 1))
fit1 <- lm(ys ~ x)
fit2 <- lm(z ~ x)
list(bys = fit1$coefficients["x",], bz = fit2$coefficients["x"])
}
fits <- do.call(rbind, lapply(1:50, function(x) sim3(bs_notfitted, corr_mat_ys)))
bys <- do.call(rbind, fits[,"bys"])
bzs <- do.call(rbind, fits[,"bz"])
coef_fit <- rlm(bzs ~ bys, maxit = 100)
par(mar = c(4,6,4,4))
plot(bs_notfitted, y = coef_fit$coefficients[-1], cex.lab = 1.5, pch = 19, col = adjustcolor(1, 0.5), cex = 1.5,
xlab = latex2exp::TeX("true value $\\beta_i$"), ylab = latex2exp::TeX("estimated value $\\hat{\\beta_i}$ from z ~ $\\Sigma y_i$"))
for(i in 1:length(coef_fit$coefficients[-1])){
segments(x0 = bs_notfitted[i], y0 = coef_fit$coefficients[-1][i] + 2*summary(coef_fit)$coefficients[,"Std. Error"][-1][i],
x1 = bs_notfitted[i], y1 = coef_fit$coefficients[-1][i] - 2*summary(coef_fit)$coefficients[,"Std. Error"][-1][i])}
abline(0,1, col = 2, lty = 2)
legend(lwd = 1, x = "topleft", legend = c("1-to-1 line", "±2SE"), lty = c(2,1), col = c(2,1))
gives
Though some of this niceness is because the sd_ys are not consistent across each of n_y; if I fix the former to, say, 5, and also fix my x to be constant across replicates and also MVN with strong indexical autocorrelation structure (e.g. using cumsum(rnorm(n_obs, sd = sqrt(1 / n_obs))) + rev(cumsum(rnorm(n_obs, sd = sqrt(1 / n_obs))))) I get:
