# MLE with unbalanced system of regressions

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?

• combine the data sets and use dummy variables to control which parameters/predictor variables are used for each response type? – Ben Bolker Jul 10 '18 at 16:26
• @BenBolker In that case the error variances would be assumed equal, an assumption I can’t make. – Jean-Paul Jul 10 '18 at 16:35
• Ben Bolker's suggestion is easy to implement in R. You can use the glmmTMB package to specify a model for the error variance. Alternatively, the weights = option with varIdent() under gls() in nlme is an option. Also, mle() in stats4 or mle2() in bbmle can handle this if you specify a model for the standard deviation. All these methods are equivalent, say except gls() which uses REML by default. – Heteroskedastic Jim Jul 11 '18 at 12:32

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.

• Absolutely awesome answer! Thanks for not only answering the question but also providing a walk-through of how to select the best final model. I am definitely going to consider this approach, thanks a lot! – Jean-Paul Jul 16 '18 at 13:46
• @Jean-Paul nice to see you found it useful! – Heteroskedastic Jim Jul 19 '18 at 0:07