# Can I use fits from two separate linear models to say something about a third?

Hi everyone :] first time posting here. I have a question about linear models. Suppose I have two models of the form:

1. $$y_1 = \beta_1x + error_1$$

2. $$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:

1. $$y_2 = \beta_4x + error_3$$

When I regress $$\beta_4$$ ~ $$\beta_1$$, ie fit:

1. $$\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:

1. $$y_2 = \beta_2x + \beta_3(\beta_1x + error_1) + error_2$$

2. $$y_2 = \beta_2x + \beta_3\beta_1x + \beta_3error_1 + error_2$$

3. $$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: • At the end,$y_2 = \beta_2\cdot x + (\beta_1\cdot x + \delta_1) + \delta_2 = B\cdot x + \Delta$? so, you need to do a linear estimation with$x$and$y_1$, and a linear estimation of$x$and$y_2$. You will find$B, \Delta$and then you can estimate$\beta_2$from$B = \beta_2+\beta_1$and$\Delta = \delta_2 + \delta_1\cdot\$
– L F
Aug 25, 2021 at 22:15