0
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ combine the data sets and use dummy variables to control which parameters/predictor variables are used for each response type? $\endgroup$ – Ben Bolker Jul 10 '18 at 16:26
  • 1
    $\begingroup$ @BenBolker In that case the error variances would be assumed equal, an assumption I can’t make. $\endgroup$ – Jean-Paul Jul 10 '18 at 16:35
  • $\begingroup$ 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. $\endgroup$ – Heteroskedastic Jim Jul 11 '18 at 12:32
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Jean-Paul Jul 16 '18 at 13:46
  • $\begingroup$ @Jean-Paul nice to see you found it useful! $\endgroup$ – Heteroskedastic Jim Jul 19 '18 at 0:07

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.