2
$\begingroup$

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:

enter image description here

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)

enter image description here

(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))

enter image description here

(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

enter image description here

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:

enter image description here

$\endgroup$
1
  • $\begingroup$ 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$ $\endgroup$ Aug 25 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.