MLE with unbalanced system of regressions

Question

I want to estimate the following system of regressions simultaneously:

$$ \begin{align} y_1 &=\alpha_1 + \beta\ x_1 + \gamma\ z_2 + \epsilon_1 \\ y_2 &=\alpha_2 + \beta\ x_2 + \gamma\ z_1 + \epsilon_2 \\ \end{align} $$

Note that the coefficients $\beta$ and $\gamma$ are the same in both regressions. Everything else is different.

My problem: the number of observations are different per regression (let's say $n_1$ and $n_2$), so how can I estimate this system of regressions using MLE?

My current solution: maximize the sum of the two individual log likelihoods

E.g.: $$ \begin{align} \text{minimize} \left[ -\log L_1 - \log L_2 \right] &= \left( \frac{n_1}{2} \log 2 \pi + \frac{n_1}{2} \log \sigma_1^2 + \frac{1}{2\sigma_1^2}\sum_{i=1}^{n_1}(y_1 - \widehat{y}_1)^2 \right) \\ &\quad + \left( \frac{n_2}{2} \log 2 \pi + \frac{n_2}{2} \log \sigma_2^2 + \frac{1}{2\sigma_2^2}\sum_{i=1}^{n_2}(y_2 - \widehat{y}_2)^2 \right) \end{align} $$

But I'm not sure this makes sense. Are there any other common ways of estimating this system?

---

_{Note: my actual problem contains four equations with various cross-equation ($\beta_i = \beta_j$) restrictions, but the question remains the same - how to estimate such a system when the observations (and error terms) are different per regression?}

Answer 1

You can stack both datasets. Then permit the variances to be heteroskedastic using generalized least squares. It is possible to accomplish this in R:

set.seed(500) # For reproducibility
library(nlme) # For heterosked variances

# Create data 1
n <- 250 # sample size 250
xa <- rnorm(n)
xb <- rbinom(n, 1, .5)
dat1 <- data.frame(
  xa, xb, y = 5 + .5 * xa + 1 * xb + rnorm(n, sd = 1))
rm(n, xa, xb)

# Create data 2
n <- 200
xa <- rnorm(n)
xb <- rbinom(n, 1, .5)
dat2 <- data.frame(
  xa, xb, y = .5 * xa + 1 * xb + rnorm(n, sd = 3))
rm(n, xa, xb)

# stack both datasets
dat <- rbind(dat1, dat2)
# Create identifier of source data
dat$id <- factor(c(rep(1, 250), rep(2, 200)))

# Constraining both sets of coefficients to be same
# using varIdent to permit heteroskedasticity
(m0 <- gls(
  y ~ 0 + id + xa + xb, dat, weights = varIdent(form = ~ 1 | id)))
Generalized least squares fit by REML
  Model: y ~ 0 + id + xa + xb 
  Data: dat 
  Log-restricted-likelihood: -857.7845

Coefficients:
        id1         id2          xa          xb 
 5.15612937 -0.06154803  0.42766940  0.74609548 

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | id 
 Parameter estimates:
       1        2 
1.000000 3.085376 
Degrees of freedom: 450 total; 446 residual
Residual standard error: 0.978681 

We can also estimate a model where we free the coefficients to be different. We use an interaction with the data source variable to enter variables as predictors of the response we want:

(m1 <- gls(
  y ~ 0 + id + id:xa + id:xb, dat, weights = varIdent(form = ~ 1 | id)))
Generalized least squares fit by REML
  Model: y ~ 0 + id + id:xa + id:xb 
  Data: dat 
  Log-restricted-likelihood: -856.9092

Coefficients:
       id1        id2     id1:xa     id2:xa     id1:xb     id2:xb 
 5.1733701 -0.2835926  0.4495109  0.1583834  0.7119283  1.1176506 

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | id 
 Parameter estimates:
       1        2 
1.000000 3.080887 
Degrees of freedom: 450 total; 444 residual
Residual standard error: 0.9785631 

Finally, we can compare both models:

anova(m0, m1)
   Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m0     1  6 1727.569 1752.171 -857.7845                        
m1     2  8 1729.818 1762.585 -856.9092 1 vs 2 1.750713  0.4167
Warning message:
In nlme::anova.lme(object = m0, m1) :
  fitted objects with different fixed effects. REML comparisons are not meaningful.

The model comparison suggests that the simpler model provides a good enough approximation to data compared to the more complicated model. The warning about REML comparisons may not be consequential. You can instead run:

anova(update(m0, method = "ML"), update(m1, method = "ML"))

In this example, the conclusions do not change when you use maximum likelihood instead of restricted maximum likelihood. I'd stick with the REML results for m0 as my final model.