$O(n_{u}^{3})$.

(a) Fixed $n_{u}$ and $n_{o}$

(b) Fixed $n$

Figure S1.1: Time needed of various matrix operations in the EM algorithm in
log10(time) in seconds versus log10($n$) in the left panel and log10($n_{u}$)
in the right panel

(a) Fixed $n_{u}$ and $n_{o}$

(b) Fixed $n$

Figure S1.2: Time needed of various matrix operations in the marginal ML
method in log10 (time) in seconds versus log10($n$) in the left panel and
lg10($n_{o}$) in the right panel.

Figure S1.1(a) illustrates the relationship between the log of the average
computing time needed to compute individual terms in the EM algorithm and the
log of the total number of units ($n$). Figure S1.1(b) shows the relationship
between the log of the average computing time and the log of the number of
unobserved units ($n_{u}$). The estimated complexities (empirical
complexities), discussed in Section S1.3, are derived by estimating the slopes
of these graphs. See Section 4.1 in the main paper for further details.
Similarly, Figure S1.2(a) displays the logarithm of the average computing time
of individual terms in the marginal ML method against the logarithm of $n$.
Figure S1.2(b) illustrates the relationship between the logarithm of average
computing time and the logarithm of the number of observed units, $n_{o}$. The
estimated complexities, discussed in Section S1.2, are determined by
estimating the slopes of these graphs.

###  S1.4 Complexity summary

Table S1.1: Theoretical and empirical (estimated) computational complexities of calculating $\textbf{V}_{oo}$ using the Direct method. Term | Description |  Theoretical complexity |  Empirical complexity
---|---|---|---
$\textbf{A}^{\top}\textbf{A}$ |  two $n\times n$ sparse matrix multiplication | $O(n)$ | $O(n^{1.032})$
$(\textbf{A}^{\top}\textbf{A})^{-1}$ |  Direct inversion of $n\times n$ sparse matrix | - | $O(n^{2.057})$
$\textbf{V}_{oo}=\left(\textbf{I}_{n}+\theta(\textbf{A}^{\top}\textbf{A})^{-1}\right)_{oo}$ |  No matrix operations | - | -

Table S1.1 presents a summary of the computational complexities associated
with each operation necessary for calculating $\textbf{V}_{oo}$ in the
marginal ML method using the direct method. Notably, the computational
complexity primarily relies on the inversion of $\textbf{A}^{\top}\textbf{A}$.
However, it is important to note that evaluating this complexity theoretically
is challenging. As a result, we approximate its complexity using simulations
as $O(n^{2.057})$, which can be expressed as approximately $O(n^{2})$.

Table S1.2: Theoretical and empirical (estimated) computational complexities of calculating $\textbf{V}_{oo}$ using the parameterisation method. Term | Description |  Theoretical complexity |  Empirical complexity
---|---|---|---
fixed $n_{o}$ | fixed $n$
$\textbf{A}^{-1}_{\textbf{B}_{o}}=\texttt{solve}((\textbf{A}^{\top})^{-1},\textbf{B}_{o}^{\top})$ |  solving a system of two sparse matrices | - | $O(n^{1.522})$ | $O(n_{o}^{0.973})$
$\textbf{A}\textbf{A}^{\top}$ |  Two $n\times n$ sparse matrix multiplication | $O(n)$ | $O(n^{1.032})$ | constant
$\textbf{R}=(\textbf{A}^{-1}_{\textbf{B}_{o}})^{\top}(\textbf{A}\textbf{A}^{\top}+\theta\textbf{I})$ |  $n_{o}\times n$ times $n\times n$ | $O(nn_{o})$ | $O(n^{1.298})$ | $O(n_{o}^{1.215})$
$\textbf{V}_{oo}=\textbf{R}\textbf{A}^{-1}_{\textbf{B}_{o}}$ |  $n_{o}\times n$ times $n\times n_{o}$ | $O(nn_{o}^{2})$ | $O(n^{1.045})$ | $O(n_{o}^{1.795})$

Table S1.2 provides a summary of the computational complexity associated with
each operation required to calculate $\textbf{V}_{oo}$ in the marginal ML
method using the parameterisation method. With the exception of the first
operation, we have access to theoretical complexity values for all other
operations, and these values align with the estimated complexities. For the
first operation, the computational complexity estimate, based on simulations,
is $O(n^{1.522}n_{o}^{0.973})$, which can be further simplified to
approximately $O(n^{1.5}n_{o})$. In Section 4.1.1 of the main paper, we
provide a comprehensive discussion on the overall complexity of calculating
$\textbf{V}_{oo}$, both in a general context and in the context of two
specific real-world scenarios.

Table S1.3: Theoretical and empirical (estimated) computational complexities of calculating T in the EM algorithm. Term | Description |  Theoretical complexity |  Empirical complexity
---|---|---|---
fixed $n_{u}$ | fixed $n$
$\textbf{A}_{\textbf{B}_{u}}=\textbf{A}\textbf{B}_{u}^{\top}$ |  Two sparse matrix multiplication | $O(nn_{u})$ | $O(n^{0.7831})$ |
$\textbf{A}\textbf{A}^{\top}$ |  Two $n\times n$ sparse matrix multiplication | $O(n)$ | $O(n^{1.032})$ | constant
$\textbf{L}_{\textbf{A}\textbf{A}^{\top}}$ |  Cholesky of local neighbourhood | $O(n^{1.5})$ |  | constant
$\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}}$ | Cholesky update | - |  | constant
$\textbf{S}=\texttt{solve}(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}},\textbf{A}_{\textbf{B}_{u}})$ |  solving a system of two sparse matrices using Cholesky factors | - | $O(n^{1.497})$ | $O(n_{u}^{0.9098})$
${\textbf{A}_{\textbf{B}_{o}}}^{\top}=\textbf{B}_{o}\textbf{A}^{\top}$ |  $n_{o}\times n$ times $n\times n$ | $O(nn_{o})$ | $O(n^{0.9477})$ |
$\textbf{M}_{uu}={\textbf{A}_{\textbf{B}_{u}}}^{\top}\textbf{S}$ |  $n_{u}\times n$ times $n\times n_{u}$ | $O(n_{u}^{2})$ | constant | $O(n_{u}^{2.1683})$
$\textbf{M}_{ou}={\textbf{A}_{\textbf{B}_{o}}}^{\top}\textbf{S}$ |  $n_{o}\times n$ times $n\times n_{u}$ | $O(n_{u}n_{o})$ | constant |
T= solve$(\textbf{M}_{uu},\textbf{M}_{uo}(\textbf{z}_{o}-\boldsymbol{\mu}_{o})$) =$\textbf{M}_{uu}^{-1}\textbf{M}_{uo}(\textbf{z}_{o}-\boldsymbol{\mu}_{o})$ |  solving a system with a dense matrix and a vector | $O(n_{u}^{3})$ | constant |

Similarly, Table S1.3 summarises the computational complexities attached to
individual operations in calculating the conditional mean
$\boldsymbol{\mu}_{u\mid o}$ for the proposed EM algorithm. Except for the
fourth and fifth operations (calculating the Cholesky update and S), we have
theoretical complexity values. We can approximate the complexity of S as
$O(n^{1.497}n_{u}^{0.9098})$, which can be expressed as approximately
$O(n^{1.5}n_{o})$. Moreover, when performing the Cholesky update operation, we
rely on the complexity findings presented by Suesse, 2018a . These findings
suggest an approximate computational complexity of $O(n^{1.5})$. A
comprehensive overview of the overall computational complexity involved in
calculating $\boldsymbol{\mu}_{u\mid o}$ is available in Section 4.1.2 of the
main paper.

##  S2 Extended simulation results

Tables S2.1 to S2.12 present the empirical means, mean squared errors of
parameter estimates and their coverages obtained from $2500$ simulated
datasets (each with $n=625$ units) for the H-SEM and H-SAM. The synthetic data
for these simulations are generated under settings similar to those employed
in the simulation study outlined in Section 4.2 of the main paper. The
coverage calculation is determined using $95\%$ Wald-type confidence
intervals, employing the formula $\hat{\delta}\pm
1.96\sqrt{\widehat{\text{Var}}(\hat{\delta})}$, where
$\delta\in\boldsymbol{\phi}$ and $\widehat{\text{Var}}(\hat{\delta})$ is the
estimated variance of the parameter $\delta$. The estimated variances for the
parameters are computed based on the formula given in Section S6 of the online
supplement.

Table S2.1: H-SEM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 10% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
$\beta_{0}$ | 1 | 0.9929 | 0.0415 | 0.9329 | 0.9932 | 0.0411 | 0.9399 | 0.9932 | 0.0411 | 0.9399
$\beta_{1}$ | 5 | 4.9962 | 0.0050 | 0.9717 | 4.9961 | 0.0050 | 0.9717 | 4.9961 | 0.0050 | 0.9717
$\rho$ | 0.8 | 0.7631 | 0.0080 | 0.9399 | 0.7779 | 0.0060 | 0.9470 | 0.7779 | 0.0060 | 0.9505
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 2.0128 | 0.1807 | 0.9364 | 1.8777 | 0.1995 | 0.9682 | 1.8777 | 0.1994 | 0.9717
$\sigma^{2}_{\textbf{y}}$ | 1 | 1.0337 | 0.2190 | 0.9011 | 1.1342 | 0.2360 | 0.9399 | 1.1342 | 0.2359 | 0.9435

For H-SEM, when estimating parameters with a small percentage of missing data,
specifically $10\%,20\%,$ and $30\%$ (see Tables S2.1, S2.3, and S2.5), we
observe that, the mean estimates obtained from all three approaches (the ODM
specification, marginal ML method, and EM algorithm) are relatively close to
the true values. However, the mean squared errors (MSE) of estimates from the
ODM specification tend to be generally higher compared to those from the CLM.

Table S2.2: H-SAM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 10% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 1.2403 | 0.0746 | 0.3307 | 0.9945 | 0.0049 | 0.9402 | 0.9945 | 0.0049 | 0.9402
${\beta}_{1}$ | 5 | 5.0602 | 0.0138 | 0.8924 | 4.9970 | 0.0070 | 0.9323 | 4.9970 | 0.0070 | 0.9323
${\rho}$ | 0.8 | 0.7519 | 0.0026 | 0.0737 | 0.8009 | 0.0001 | 0.9422 | 0.8009 | 0.0001 | 0.9422
${\sigma}^{2}_{\boldsymbol{\epsilon}}$ | 2 | 1.8826 | 0.1091 | 0.9402 | 2.0182 | 0.0686 | 0.9502 | 2.0182 | 0.0686 | 0.9502
$\sigma^{2}_{\textbf{y}}$ | 1 | 2.1226 | 1.3899 | 0.0319 | 0.9621 | 0.0385 | 0.9183 | 0.9621 | 0.0384 | 0.9183

For H-SAM, when dealing with $10\%$, $20\%$, $30\%$, $70\%$, and $80\%$
missing data, the marginal ML method and EM algorithm consistently provide
nearly identical estimates with a coverage probability of approximately 0.95
(i.e., $95\%$) and exhibit low mean squared errors (MSE). In contrast, the
estimates obtained from the observed data model specification deviate
significantly from the true values, displaying higher MSE compared to the
marginal ML method and EM algorithm. Some parameters even exhibit 0 coverage
in certain cases. See Tables S2.2, S2.4, S2.6, S2.8, and S2.10 for further
details.

Table S2.3: H-SEM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 20% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
$\beta_{0}$ | 1 | 1.0422 | 0.0494 | 0.8774 | 1.0429 | 0.0489 | 0.9151 | 1.0429 | 0.0489 | 0.9151
$\beta_{1}$ | 5 | 5.0062 | 0.0066 | 0.9623 | 5.0060 | 0.0066 | 0.9623 | 5.0060 | 0.0066 | 0.9623
$\rho$ | 0.8 | 0.7400 | 0.0138 | 0.9214 | 0.7743 | 0.0077 | 0.9308 | 0.7743 | 0.0077 | 0.9245
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 2.1273 | 0.2788 | 0.8679 | 1.8771 | 0.2538 | 0.9434 | 1.8770 | 0.2537 | 0.9403
$\sigma^{2}_{\textbf{y}}$ | 1 | 0.9719 | 0.2932 | 0.8270 | 1.1424 | 0.2631 | 0.9308 | 1.1424 | 0.2630 | 0.9245
Table S2.4: H-SAM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 20% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 1.5631 | 0.3570 | 0.0351 | 1.0041 | 0.0045 | 0.9396 | 1.0041 | 0.0045 | 0.9396
${\beta}_{1}$ | 5 | 5.1861 | 0.0489 | 0.6316 | 5.0077 | 0.0071 | 0.9688 | 5.0077 | 0.0071 | 0.9688
$\rho$ | 0.8 | 0.6892 | 0.0129 | 0.0000 | 0.7994 | 0.0001 | 0.9376 | 0.7994 | 0.0001 | 0.9376
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 1.8075 | 0.2354 | 0.9240 | 2.0213 | 0.0833 | 0.9532 | 2.0213 | 0.0832 | 0.9532
$\sigma^{2}_{\textbf{y}}$ | 1 | 3.4897 | 6.5525 | 0.0000 | 0.9658 | 0.0454 | 0.9181 | 0.9659 | 0.0454 | 0.9181
Table S2.5: H-SEM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 30% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 0.9690 | 0.0457 | 0.8806 | 0.9659 | 0.0454 | 0.9463 | 0.9659 | 0.0454 | 0.9463
${\beta}_{1}$ | 5 | 4.9927 | 0.0083 | 0.9522 | 4.9926 | 0.0081 | 0.9493 | 4.9926 | 0.0081 | 0.9493
$\rho_{0}$ | 0.8 | 0.7139 | 0.0182 | 0.9373 | 0.7817 | 0.0067 | 0.9612 | 0.7817 | 0.0067 | 0.9582
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 2.2285 | 0.2808 | 0.8627 | 1.8757 | 0.2600 | 0.9612 | 1.8757 | 0.2598 | 0.9582
$\sigma^{2}_{\textbf{y}}$ | 1 | 0.9166 | 0.2770 | 0.8507 | 1.1203 | 0.2700 | 0.9403 | 1.1203 | 0.2699 | 0.9403
Table S2.6: H-SAM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 30% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 1.9794 | 1.0589 | 0.0054 | 1.0039 | 0.0046 | 0.9458 | 1.0039 | 0.0046 | 0.9458
${\beta}_{1}$ | 5 | 5.3251 | 0.1262 | 0.3469 | 4.9994 | 0.0079 | 0.9485 | 4.9994 | 0.0079 | 0.9485
$\rho$ | 0.8 | 0.6136 | 0.0360 | 0.0000 | 0.7998 | 0.0001 | 0.9512 | 0.7998 | 0.0001 | 0.9512
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 1.5454 | 0.6575 | 0.9079 | 2.0069 | 0.1144 | 0.9377 | 2.0069 | 0.1142 | 0.9377
$\sigma^{2}_{\textbf{y}}$ | 1 | 5.3717 | 20.1868 | 0.0000 | 0.9734 | 0.0495 | 0.9350 | 0.9734 | 0.0494 | 0.9350
Table S2.7: H-SEM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 70% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 0.9863 | 0.0809 | 0.7903 | 0.9797 | 0.0820 | 0.9032 | 0.9796 | 0.0820 | 0.9032
${\beta}_{1}$ | 5 | 4.9916 | 0.0147 | 0.9886 | 4.9897 | 0.0153 | 0.9839 | 4.9897 | 0.0153 | 0.9839
$\rho$ | 0.8 | 0.5302 | 0.1331 | 0.6613 | 0.7726 | 0.0136 | 0.9032 | 0.7726 | 0.0135 | 0.9032
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 2.0334 | 1.8126 | 0.8065 | 1.7241 | 1.0984 | 0.9677 | 1.7248 | 1.0915 | 0.9677
$\sigma^{2}_{\textbf{y}}$ | 1 | 1.6499 | 2.4837 | 0.8065 | 1.2912 | 0.8625 | 0.8548 | 1.2901 | 0.8574 | 0.8548

In the case of fitting the H-SEM with $70\%$ missing values, the MSEs are
lower, and the coverage is generally higher for estimates derived from the CLM
specification compared to estimates from the ODM specification. See Table S2.7
for further details. In contrast to scenarios with lower missing value
percentages such as $10\%$, $20\%$, and $30\%$, the estimate for $\rho$ from
ODM considerably deviates from the true value.

Table S2.8: H-SAM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 70% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 3.6927 | 7.5528 | 0.0000 | 1.0044 | 0.0073 | 0.9386 | 1.0044 | 0.0073 | 0.9386
${\beta}_{1}$ | 5 | 5.9984 | 1.0796 | 0.0702 | 5.0053 | 0.0212 | 0.9474 | 5.0053 | 0.0212 | 0.9474
$\rho$ | 0.8 | 0.3444 | 0.2095 | 0.0000 | 0.7987 | 0.0002 | 0.9561 | 0.7987 | 0.0002 | 0.9561
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 0.0387 | 3.8995 | 0.9754 | 2.0980 | 0.4133 | 0.9649 | 2.0970 | 0.4096 | 0.9649
$\sigma^{2}_{\textbf{y}}$ | 1 | 15.1225 | 205.5219 | 0.0526 | 0.9440 | 0.1491 | 0.9123 | 0.9444 | 0.1482 | 0.9123
Table S2.9: H-SEM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 80% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 0.9424 | 0.0583 | 0.9000 | 0.9522 | 0.0514 | 0.9734 | 0.9522 | 0.0514 | 0.9734
${\beta}_{1}$ | 5 | 4.9346 | 0.0306 | 0.9333 | 4.9342 | 0.0323 | 0.9333 | 4.9343 | 0.0323 | 0.9333
$\rho$ | 0.8 | 0.5012 | 0.1742 | 0.8000 | 0.8162 | 0.0105 | 0.8333 | 0.8159 | 0.0105 | 0.8333
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 2.1295 | 1.3991 | 0.8333 | 1.9895 | 0.8954 | 0.8667 | 1.9879 | 0.8903 | 0.8667
$\sigma^{2}_{\textbf{y}}$ | 1 | 1.3452 | 2.1345 | 0.8000 | 0.8016 | 0.4482 | 0.8000 | 0.8021 | 0.4469 | 0.8000

Fitting the H-SEM with $80\%$ missing values using the ODM specification leads
to inaccurate estimates compared to the CLM specification. The ODM
specification has higher MSE values compared to those of CLM for the
parameters $\rho$, ${\sigma}^{2}_{\boldsymbol{\epsilon}}$ and
$\sigma^{2}_{\textbf{y}}$, as displayed in Table S2.9.

Table S2.10: H-SAM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 80% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 4.0162 | 9.5202 | 0.0000 | 0.9838 | 0.0117 | 0.8605 | 0.9838 | 0.0117 | 0.8605
${\beta}_{1}$ | 5 | 6.0322 | 1.2216 | 0.2093 | 4.9542 | 0.0420 | 0.9302 | 4.9543 | 0.0420 | 0.9302
$\rho$ | 0.8 | 0.2737 | 0.2799 | 0.0000 | 0.8017 | 0.0003 | 0.9070 | 0.8017 | 0.0003 | 0.9070
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 0.3197 | 3.7892 | 0.9877 | 2.2445 | 0.7745 | 0.9535 | 2.2406 | 0.7614 | 0.9535
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 1 | 16.2810 | 241.6783 | 0.6047 | 0.8367 | 0.2287 | 0.8837 | 0.8384 | 0.2257 | 0.9070

In Section 3.3 of the main paper, we discuss the limitations associated with
employing the EM algorithm when dealing with large values of $n_{u}$. When
dealing with datasets with $90\%$ missing values, the EM algorithm exhibits
slow convergence. As a result, we only compare the estimates from the ODM with
the estimates from the marginal ML method.

Tables S2.11 and S2.12 present the mean estimates, the MSE, and the coverage
of estimates obtained from the ODM and CLM estimated using the marginal ML
method with 90% missing values.

Table S2.11: H-SEM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 90% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Entries marked with ’NC’ indicate that the EM algorithm failed to converge for a significant number of simulations, preventing the reliable calculation of mean estimates, MSE, and coverages. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 0.9966 | 0.1101 | 0.8848 | 0.9939 | 0.1077 | 0.9168 | NC | NC | NC
${\beta}_{1}$ | 5 | 5.0103 | 0.0738 | 0.8542 | 5.0141 | 0.0706 | 0.9424 | NC | NC | NC
${\rho}$ | 0.8 | 0.2781 | 0.4004 | 0.6112 | 0.5547 | 0.2621 | 0.8384 | NC | NC | NC
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 1.1193 | 3.2489 | 0.9360 | 0.9749 | 2.5725 | 0.9808 | NC | NC | NC
$\sigma^{2}_{\textbf{y}}$ | 1 | 2.7998 | 6.2493 | 0.8588 | 1.9075 | 2.4003 | 0.9040 | NC | NC | NC

The MSEs of estimates obtained from the CLM specifications are consistently
lower than those from the ODM specification for the H-SEM with $90\%$ missing
values, as illustrated in Table S2.11. These results imply that the CLM
specification displays greater robustness to high percentages of missing data,
resulting in more accurate estimates for the H-SEM.

Table S2.12: H-SAM mean estimates (est.), mean squared error (m.s.e.), and coverage (cov.) of estimated parameters by the ODM, the marginal ML method, and the EM algorithm - 90% of missing data. These attributes are computed from 2500 simulated datasets, each with $n=625$ units. Entries marked with ’NC’ indicate that the EM algorithm failed to converge for a significant number of simulations, preventing the reliable calculation of mean estimates, MSE, and coverages. Parameter |  True value | ODM | CLM
---|---|---|---
ML | Marginal ML | EM algorithm
|  | est. | m.s.e | cov. | est. | m.s.e | cov. | est. | m.s.e | cov.
${\beta}_{0}$ | 1 | 4.4343 | 12.6305 | 0.0045 | 1.0019 | 0.0209 | 0.9189 | NC | NC | NC
${\beta}_{1}$ | 5 | 6.2834 | 2.0270 | 0.4009 | 5.0123 | 0.0812 | 0.9279 | NC | NC | NC
$\rho$ | 0.8 | 0.2849 | 0.2750 | 0.0180 | 0.7991 | 0.0007 | 0.9009 | NC | NC | NC
$\sigma^{2}_{\boldsymbol{\epsilon}}$ | 2 | 3.7648 | 47.2668 | 0.9775 | 2.2564 | 1.5172 | 0.9910 | NC | NC | NC
$\sigma^{2}_{\textbf{y}}$ | 1 | 13.8161 | 213.6280 | 0.9640 | 0.7381 | 0.4297 | 0.9595 | NC | NC | NC

Similar to the scenarios with missing percentages of $10\%$, $20\%$, $30\%$,
$70\%$, and $80\%$, the CLM specification consistently yields estimates that
are closer to the true parameters, exhibiting lower MSEs and higher coverages
compared to ODM for H-SAM with $90\%$ of missing values.

In summary, we find that estimating the H-SAM with ODM specification results
in inaccurate parameter estimates for any level of missing value percentage,
whereas CLM consistently produces highly accurate parameter estimates.
Consequently, we recommend adopting the CLM specification for estimating H-SAM
regardless of the missing value percentage. On the other hand, for H-SEM, when
dealing with a low percentage of missing values ($10\%,20\%,30\%$), both ODM
and CLM provide accurate parameter estimates. However, as the missing value
percentage increases to $70\%,80\%$, and $90\%$, ODM produces estimates that
are considerably different from the true values. In contrast, CLM yields
estimates closer to the true values, demonstrating lower MSE, particularly for
parameters such as $\rho$, ${\sigma}^{2}_{\boldsymbol{\epsilon}}$ and
$\sigma^{2}_{\textbf{y}}$. Therefore, we recommend utilising CLM for H-SEM,
especially in scenarios with higher percentages of missing values.

##  S3 Extended Real data analysis results

Table S3.1 displays estimates and standard errors for fixed effects
($\boldsymbol{\beta}$’s), $\rho$, $\sigma^{2}_{\epsilon}$, and
$\sigma^{2}_{\textbf{y}}$ for the H-SAM with $10\%$ missing values obtained
using FDM, ODM, and CLM estimated using the marginal ML method and the EM
algorithm. The table includes the marginal log-likelihood values for the
marginal ML and the EM methods, the computing time, the time taken to compute
the standard errors of the parameter estimates, and the number of iterations
for the EM algorithm.

According to the results presented in the complexity analysis in Section 4.1
of the main paper, when dealing with low percentages of missing values, the EM
algorithm is generally expected to offer computational advantages over the
marginal ML method for both H-SEM and H-SAM. However, contrary to this
expectation, as indicated in the third last row of Table S3.1, the EM
algorithm takes longer to converge than the marginal ML method. This
observation can be attributed to the high number of iterations (10) in the EM
algorithm. It is important to note that despite the difference in computation
time, both marginal ML and EM algorithm methods produce identical parameter
estimates. Furthermore, parameter estimates obtained from CLM are found to be
closer to the estimates obtained from the FDM than those from the ODM, which
relies solely on observed data.

Table S3.4 displays estimates and standard errors for fixed effects
($\boldsymbol{\beta}$’s), $\rho$, $\sigma^{2}_{\epsilon}$, and
$\sigma^{2}_{\textbf{y}}$ for the H-SEM with $90\%$ missing values obtained
using FDM, ODM, and CLM estimated using the marginal ML method and the EM
algorithm. The table includes the marginal log-likelihood values for the
marginal ML and the EM methods, the computing time, the time taken to compute
the standard errors of the parameter estimates, and the number of iterations
for the EM algorithm.

It is important to highlight that, similar to the case of fitting H-SAM with
$90\%$ missing data as discussed in Section 5 of the main paper, we terminated
the EM algorithm for H-SEM with $90\%$ missing data after 100 iterations
(after about 49 hours and 53 minutes) due to computational constraints. At
this point, the EM algorithm had not reached convergence, and its execution
time was significantly longer in comparison to the marginal ML method. This
observation aligns with the conclusions drawn from our computational
complexity study in Section 4.1 of the main paper, emphasising that when faced
with a substantial proportion of missing data, the marginal ML method
demonstrates superior computational efficiency compared to the EM algorithm.
Note that, in this example, the computational cost of the marginal ML and EM
algorithms is large because the total number of observations (n) is large.

Table S3.1: Parameter estimates (est.) with their standard errors (S.E.), marginal log-likelihood values for the marginal ML and the EM algorithms, estimation time and the time to calculate the standard errors of the parameters (Time std.error) in seconds for fitting H-SAM using FDM-ML for the full data set, and ODM-ML, CLM-Marginal ML, and CLM-EM for a partially observed data set with $10\%$ of missing data, and the number of EM iterations. | FDM-ML | ODM-ML | CLM-Marginal ML | CLM-EM
---|---|---|---|---
est. | S.E. | est. | S.E. | est. | S.E. | est. | S.E.
Intercept | -0.1124 | 0.0507 | 2.7726 | 0.0895 | -0.0962 | 0.0538 | -0.0962 | 0.0538
age | 0.9565 | 0.0429 | 1.9633 | 0.0805 | 0.9512 | 0.0453 | 0.9512 | 0.0453
$\text{age}^{2}$ | -1.5790 | 0.0797 | -4.0116 | 0.1458 | -1.5702 | 0.0848 | -1.5701 | 0.0848
$\text{age}^{3}$ | 0.3697 | 0.0440 | 1.2107 | 0.0789 | 0.3714 | 0.0458 | 0.3714 | 0.0458
log(lotsize) | 0.0413 | 0.0022 | 0.1724 | 0.0043 | 0.0400 | 0.0024 | 0.0400 | 0.0024
rooms | -0.0052 | 0.0026 | 0.0090 | 0.0044 | -0.0045 | 0.0028 | -0.0045 | 0.0028
log(TLA) | 0.4454 | 0.0083 | 0.8930 | 0.0141 | 0.4351 | 0.0098 | 0.4351 | 0.0098
beds | 0.0129 | 0.0039 | -0.0174 | 0.0065 | 0.0111 | 0.0041 | 0.0111 | 0.0041
syear1994 | 0.0357 | 0.0066 | 0.0461 | 0.0106 | 0.0333 | 0.0069 | 0.0333 | 0.0069
syear1995 | 0.0710 | 0.0064 | 0.0835 | 0.0103 | 0.0670 | 0.0067 | 0.0670 | 0.0067
syear1996 | 0.0864 | 0.0063 | 0.1083 | 0.0100 | 0.0810 | 0.0065 | 0.0810 | 0.0065
syear1997 | 0.1191 | 0.0062 | 0.1451 | 0.0099 | 0.1134 | 0.0065 | 0.1134 | 0.0065
syear1998 | 0.1675 | 0.0064 | 0.2045 | 0.0102 | 0.1635 | 0.0067 | 0.1635 | 0.0067
$\rho$ | 0.6727 | 0.0001 | 0.0220 | 0.0300 | 0.6794 | 0.0051 | 0.6794 | 0.0051
$\sigma^{2}_{\textbf{y}}$ | 0.0399 | 0.0008 | 0.1745 | 0.0062 | 0.0399 | 0.0013 | 0.0399 | 0.0013
$\sigma^{2}_{\epsilon}$ | 0.042 | 0.0009 |  0.0001 | 0.0016 | 0.0418 | 0.0012 | 0.0418 | 0.0012
number of EM iterations |  |  |  | 10
estimation time (s) | 6.319 | 20.091 | 642.4235 | 1438.0863
Time std.error (s) | 5.926 | 9.0800 | 93.8317 | 85.8197
marginal log-likelihood |  |  | -6838.894 | -6838.894
Table S3.2:  Parameter estimates (est.) with their standard errors (S.E.), marginal log-likelihood values for the marginal ML and the EM algorithms, estimation time and the time to calculate the standard errors of the parameters (Time std.error) in seconds for fitting H-SAM using FDM-ML for the full data set, and ODM-ML, CLM-Marginal ML, and CLM-EM for a partially observed data set with $90\%$ of missing data, and the number of EM iterations. Entries labeled as ’NC’ indicate instances where the EM algorithm did not successfully converge. | FDM-ML | ODM-ML | CLM-Marginal ML | CLM-EM
---|---|---|---|---
est. | S.E. | est. | S.E. | est. | S.E. | est. | S.E.
Intercept | -0.1124 | 0.0507 | 3.2539 | 0.2683 | 0.1863 | 0.1511 | NC
age | 0.9565 | 0.0429 | 1.6800 | 0.2551 | 0.8215 | 0.1219 | NC
$\text{age}^{2}$ | -1.5790 | 0.0797 | -3.7745 | 0.4711 | -1.2960 | 0.2327 | NC
$\text{age}^{3}$ | 0.3697 | 0.0440 | 1.1916 | 0.2583 | 0.0989 | 0.1264 | NC
log(lotsize) | 0.0413 | 0.0022 | 0.1701 | 0.0124 | 0.0414 | 0.0064 | NC
rooms | -0.0052 | 0.0026 | 0.0296 | 0.0134 | -0.0042 | 0.0082 | NC
log(TLA) | 0.4454 | 0.0083 | 0.8581 | 0.0422 | 0.5234 | 0.0324 | NC
beds | 0.0129 | 0.0039 | -0.0312 | 0.0194 | -0.0143 | 0.0125 | NC
syear1994 | 0.0357 | 0.0066 | 0.0620 | 0.0311 | 0.0631 | 0.0213 | NC
syear1995 | 0.0710 | 0.0064 | 0.0900 | 0.0305 | 0.0816 | 0.0207 | NC
syear1996 | 0.0864 | 0.0063 | 0.1005 | 0.0297 | 0.1010 | 0.0201 | NC
syear1997 | 0.1191 | 0.0062 | 0.1341 | 0.0290 | 0.1288 | 0.0200 | NC
syear1998 | 0.1675 | 0.0064 | 0.1629 | 0.0304 | 0.1518 | 0.0212 | NC
$\rho$ | 0.6727 | 0.0001 | 0.0027 | 0.0311 | 0.6046 | 0.0201 | NC
$\sigma^{2}_{\textbf{y}}$ | 0.0399 | 0.0008 | 0.0882 | 0.0114 | 0.0682 | 0.0070 | NC
$\sigma^{2}_{\epsilon}$ | 0.042 | 0.0009 | 0.0882 | 0.0234 | 0.0203 | 0.0080 | NC
number of EM iterations |  |  |  | 100
estimation time (s) | 6.319 | 0.216 | 18.7909 | 101916.2943
Time std.error (s) | 5.926 | 0.836 | 46.3919 | 41.6773
marginal log-likelihood |  |  | -1120.645 | NC
Table S3.3: Parameter estimates (est.) with their standard errors (S.E.), marginal log-likelihood values for the marginal ML and the EM algorithms, estimation time and the time to calculate the standard errors of the parameters (Time std.error) in seconds for fitting H-SEM using FDM-ML for the full data set, and ODM-ML, CLM-Marginal ML, and CLM-EM for a partially observed data set with $10\%$ of missing data, and the number of EM iterations. | FDM-ML | ODM-ML | CLM-Marginal ML | CLM-EM
---|---|---|---|---
est. | S.E. | est. | S.E. | est. | S.E. | est. | S.E.
Intercept | 5.2578 | 0.0748 | 4.3654 | 0.0831 | 5.2025 | 0.0793 | 5.2024 | 0.079
age | 0.6994 | 0.0793 | 1.2398 | 0.0848 | 0.8220 | 0.0837 | 0.8203 | 0.0837
$\text{age}^{2}$ | -1.7558 | 0.1321 | -2.9150 | 0.1440 | -1.9681 | 0.1398 | -1.9681 | 0.1398
$\text{age}^{3}$ | 0.6355 | 0.0659 | 1.0907 | 0.0733 | 0.7360 | 0.0700 | 0.7348 | 0.0700
log(lotsize) | 0.1458 | 0.0046 | 0.2008 | 0.0048 | 0.1482 | 0.0049 | 0.1481 | 0.0049
rooms | 0.0056 | 0.0029 | 0.0069 | 0.0034 | 0.0062 | 0.0031 | 0.0062 | 0.0031
log(TLA) | 0.6038 | 0.0103 | 0.6606 | 0.0119 | 0.6059 | 0.0109 | 0.6059 | 0.0109
beds | 0.0164 | 0.0043 | 0.0131 | 0.0050 | 0.0164 | 0.0045 | 0.0164 | 0.0045
syear1994 | 0.0365 | 0.0067 | 0.0400 | 0.0079 | 0.0339 | 0.0072 | 0.0339 | 0.0072
syear1995 | 0.0799 | 0.0066 | 0.0805 | 0.0077 | 0.0781 | 0.0070 | 0.0781 | 0.0070
syear1996 | 0.0962 | 0.0064 | 0.1006 | 0.0075 | 0.0940 | 0.0068 | 0.0940 | 0.0068
syear1997 | 0.1413 | 0.0063 | 0.1438 | 0.0074 | 0.1393 | 0.0067 | 0.1393 | 0.0068
syear1998 | 0.1937 | 0.0065 | 0.1954 | 0.0076 | 0.1937 | 0.0069 | 0.1937 | 0.0069
$\rho$ | 0.9866 | 0.0002 | 0.8654 | 0.0049 | 0.9868 | 0.0002 | 0.9868 | 0.0002
$\sigma^{2}_{\textbf{y}}$ | 0.0004 | 0.0001 | 0.0226 | 0.0011 | 0.0004 | 0.0001 | 0.0004 | 0.0001
$\sigma^{2}_{\epsilon}$ | 0.0685 | 0.0007 | 0.0603 | 0.0012 | 0.0691 | 0.0007 | 0.0691 | 0.0007
number of EM iterations |  |  |  | 8
estimation time (s) | 14.249 | 4.963 | 1803.9157 | 1546.3783
Time std.error (s) | 5.849 | 5.52 | 59.2792 | 65.7387
marginal log-likelihood |  |  | -5885.136 | -5885.098
Table S3.4: Parameter estimates (est.) with their standard errors (S.E.), marginal log-likelihood values for the marginal ML and the EM algorithms, estimation time and the time to calculate the standard errors of the parameters (Time std.error) in seconds for fitting H-SEM using FDM-ML for the full data set, and ODM-ML, CLM-Marginal ML, and CLM-EM for a partially observed data set with $90\%$ of missing data, and the number of EM iterations. Entries labeled as ’NC’ indicate instances where the EM algorithm did not successfully converge. | FDM-ML | ODM-ML | CLM-Marginal ML | CLM-EM
---|---|---|---|---
est. | S.E. | est. | S.E. | est. | S.E. | est. | S.E.
Intercept | 5.2578 | 0.0748 | 3.3498 | 0.2640 | 3.4952 | 0.2446 | NC
age | 0.6994 | 0.0793 | 1.5135 | 0.2553 | 1.1073 | 0.2472 | NC
$\text{age}^{2}$ | -1.7558 | 0.1321 | -3.4992 | 0.4678 | -2.6208 | 0.4332 | NC
$\text{age}^{3}$ | 0.6355 | 0.0659 | 1.1014 | 0.2548 | 0.7277 | 0.2287 | NC
log(lotsize) | 0.1458 | 0.0046 | 0.1806 | 0.0125 | 0.1743 | 0.0133 | NC
rooms | 0.0056 | 0.0029 | 0.0264 | 0.0128 | 0.0021 | 0.0110 | NC
log(TLA) | 0.6038 | 0.0103 | 0.8340 | 0.0411 | 0.8300 | 0.0365 | NC
beds | 0.0164 | 0.0043 | -0.0241 | 0.0187 | -0.0049 | 0.0165 | NC
syear1994 | 0.0365 | 0.0067 | 0.0660 | 0.0299 | 0.0770 | 0.0258 | NC
syear1995 | 0.0799 | 0.0066 | 0.0862 | 0.0294 | 0.1206 | 0.0253 | NC
syear1996 | 0.0962 | 0.0064 | 0.0977 | 0.0287 | 0.1384 | 0.0247 | NC
syear1997 | 0.1413 | 0.0063 | 0.1386 | 0.0277 | 0.1744 | 0.0242 | NC
syear1998 | 0.1937 | 0.0065 | 0.1575 | 0.0293 | 0.2148 | 0.0253 | NC
$\rho$ | 0.9866 | 0.0002 | 0.6866 | 0.0012 | 0.9936 | 0.0006 | NC
$\sigma^{2}_{\textbf{y}}$ | 0.0004 | 0.0001 | 0.1643 | 0.0030 | 0.0001 | 0.0002 | NC
$\sigma^{2}_{\epsilon}$ | 0.0685 | 0.0007 | 0.0001 | 0.0081 | 0.0837 | 0.0035 | NC
number of EM iterations |  |  |  | 100
estimation time (s) | 14.249 | 0.968 | 52.02428 | 179518.2548
Time std.error (s) | 5.849 | 0.363 | 42.44297 | 44.4597
marginal log-likelihood |  |  | -1159.595 | NC

##  S4 Additional proofs for the marginal ML method

This section derives the analytical forms of ML estimators for
$\boldsymbol{\beta}$ and $\omega$ for the proposed marginal ML method
presented in Section 3.2 of the main paper. By differentiating the marginal
log-likelihood for $\textbf{z}_{o}$ in Equation (3.3) in the main paper with
respect to $\boldsymbol{\beta}$ we get

$\begin{split}\diffp{\textrm{log}f(\textbf{z}_{o};\omega,\boldsymbol{\theta},\rho,\boldsymbol{\beta})}{\boldsymbol{\beta}}&=-\frac{1}{2\omega}\diffp{(\textbf{r}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{r}_{o})}{\boldsymbol{\beta}},\end{split}$
(S4.1)

by substituting
$\textbf{r}_{o}=\textbf{z}_{o}-\tilde{\textbf{X}}_{o}\boldsymbol{\beta}$,

$\begin{split}\diffp{\textrm{log}f(\textbf{z}_{o};\omega,\boldsymbol{\theta},\rho,\boldsymbol{\beta})}{\boldsymbol{\beta}}&=-\frac{1}{2\omega}\diffp{(\textbf{z}_{o}-\tilde{\textbf{X}}_{o}\boldsymbol{\beta})^{\top}\textbf{V}_{oo}^{-1}(\textbf{z}_{o}-\tilde{\textbf{X}}_{o}\boldsymbol{\beta})}{\boldsymbol{\beta}}\\\
&=-\frac{2}{2\omega}(-\tilde{\textbf{X}}_{o})^{\top}\textbf{V}_{oo}^{-1}(\textbf{z}_{o}-\tilde{\textbf{X}}_{o}\boldsymbol{\beta})\\\
&=\frac{1}{\omega}(\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{z}_{o}-\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\tilde{\textbf{X}}_{o}\boldsymbol{\beta}).\\\
\end{split}$ (S4.2)

Then, by setting Equation (S4.2) to zero, we obtain the ML estimator for
$\boldsymbol{\beta}$

$\begin{split}0&=\frac{1}{\omega}(\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{z}_{o}-\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\tilde{\textbf{X}}_{o}\boldsymbol{\beta})\\\
\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\tilde{\textbf{X}}_{o}\boldsymbol{\beta}&=\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{z}_{o}\\\
\hat{\boldsymbol{\beta}}(\rho,\theta)&=\left(\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\tilde{\textbf{X}}_{o}\right)^{-1}\tilde{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{z}_{o}.\end{split}$
(S4.3)

Similarly, the ML estimator for $\omega$ is derived by differentiating (3.3)
with respect to $\omega$ as

$\begin{split}\diffp{\textrm{log}f(\textbf{z}_{o};\omega,\boldsymbol{\theta},\rho,\boldsymbol{\beta})}{\omega}&=-\frac{n_{o}}{2\omega}+\frac{1}{2\omega^{2}}\textbf{r}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{r}_{o}.\end{split}$
(S4.4)

By setting Equation (S4.4) to zero, the ML estimator for $\hat{\omega}$ is

$\hat{\omega}(\rho,\theta)=\frac{\textbf{r}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{r}_{o}}{n_{o}}.$
(S4.5)

##  S5 Additional proof for the EM algorithm

###  S5.1 Maximum Likelihood estimators for the parameters
$\boldsymbol{\beta}$ and ${\omega}$

This section derives the analytical forms of ML estimators for
$\boldsymbol{\beta}$ and $\omega$ for the proposed EM algorithm presented in
Section 3.3 of the main paper. By differentiating
$Q(\boldsymbol{\phi}\mid\boldsymbol{\phi}^{\prime})$ in Equation (3.15) in the
main paper with respect to $\boldsymbol{\beta}$, we obtain

$\begin{split}\diffp{Q(\boldsymbol{\phi}\mid\boldsymbol{\phi}^{\prime})}{\boldsymbol{\beta}}&=-\frac{1}{2\omega}\diffp{(\textbf{r}_{u\mid
o}^{\top}\textbf{M}\textbf{r}_{u\mid o})}{\boldsymbol{\beta}},\\\ \end{split}$
(S5.1)

and by substituting $\textbf{r}_{u\mid o}=E_{u\mid
o}(\textbf{z})-\boldsymbol{\mu}$, we obtain

$\begin{split}\diffp{Q(\boldsymbol{\phi}\mid\boldsymbol{\phi}^{\prime})}{\boldsymbol{\beta}}&=-\frac{1}{2\omega}\diffp{(E_{u\mid
o}(\textbf{z})-\boldsymbol{\mu})^{\top}\textbf{M}(E_{u\mid
o}(\textbf{z})-\boldsymbol{\mu})}{\boldsymbol{\beta}}\\\
&=-\frac{2}{2\omega}(-\tilde{\textbf{X}})^{\top}\textbf{M}^{\top}(E_{u\mid
o}(\textbf{z})-\boldsymbol{\mu})\\\
&=\frac{1}{\omega}(\tilde{\textbf{X}}^{\top}\textbf{M}E_{u\mid
o}(\textbf{z})-\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}\boldsymbol{\beta})\end{split}$
(S5.2)

Then, by setting Equation (S5.2) to zero, we obtain the ML estimator for
$\boldsymbol{\beta}$,

$\begin{split}0&=\frac{1}{\omega}(\tilde{\textbf{X}}^{\top}\textbf{M}E_{u\mid
o}(\textbf{z})-\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}\boldsymbol{\beta})\\\
\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}\boldsymbol{\beta}&=\tilde{\textbf{X}}^{\top}\textbf{M}E_{u\mid
o}(\textbf{z})\\\
\hat{\boldsymbol{\beta}}(\rho,\theta)&=\left(\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}\right)^{-1}\tilde{\textbf{X}}^{\top}\textbf{M}E_{u\mid
o}(\textbf{z}).\end{split}$ (S5.3)

Similarly, the ML estimator for $\omega$ is derived by differentiating the
conditional expectation in Equation (3.15) in the main paper with respect to
$\omega$ as

$\diffp{Q(\boldsymbol{\phi}\mid\boldsymbol{\phi}^{\prime})}{\omega}=-\frac{n}{2\omega}+\frac{1}{2\omega^{2}}\textbf{r}^{\top}_{u\mid\o}\textbf{M}\textbf{r}_{u\mid\o}.$
(S5.4)

By setting Equation (S5.4) to zero, the ML estimator for $\hat{\omega}$ is,

$\hat{\omega}(\rho,\theta)=\frac{\textbf{r}_{u\mid
o}^{\top}\textbf{M}\textbf{r}_{u\mid o}}{n}.$ (S5.5)

###  S5.2 Simplifying other terms in the EM algorithm.

This section simplifies the expressions for the terms M,
$\text{log}|\textbf{M}|$,
$\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}$, and
$\textbf{r}_{u\mid o}^{\top}\textbf{M}\textbf{r}_{u\mid o}$ found in the EM
algorithm presented in Section 3.3 of the main paper. This simplification aims
to reduce the computational burden associated with direct calculations. First,
the matrix M can be expressed as

$\begin{split}\textbf{M}=\textbf{V}^{-1}&=(\textbf{I}_{n}+\theta(\textbf{A}^{\top}\textbf{A})^{-1})^{-1}\\\
&=\textbf{A}^{\top}(\textbf{A}^{\top})^{-1}(\textbf{I}_{n}+\theta\textbf{A}^{-1}(\textbf{A}^{\top})^{-1})^{-1}\textbf{A}^{-1}\textbf{A}\\\
&=\textbf{A}^{\top}\left(\textbf{A}(\textbf{I}_{n}+\theta\textbf{A}^{-1}(\textbf{A}^{\top})^{-1})\textbf{A}^{\top}\right)^{-1}\textbf{A}\\\
&=\textbf{A}^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}.\end{split}$
(S5.6)

Second, the determinant of the $\text{log}|\textbf{M}|$ is expressed as

$\begin{split}|\textbf{M}|&=|\textbf{A}^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}|\\\
|\textbf{M}|&=|\textbf{A}^{\top}||\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}||\textbf{A}|\\\
\text{log}|\textbf{M}|&=\text{log}|\textbf{A}^{\top}|-\text{log}|\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta|+\text{log}|\textbf{A}|\\\
\text{log}|\textbf{M}|&=-\text{log}|\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}|+\text{log}|\textbf{A}\textbf{A}^{\top}|.\end{split}$
(S5.7)

We define the following two terms:

$\textbf{C}=(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\tilde{\textbf{X}}_{\textbf{A}},$
(S5.8)

and

$\textbf{c}=(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\textbf{z}_{\textbf{A}_{u\mid
o}},$ (S5.9)

where $\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}}$ is the
Cholesky factor of $\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}$,
$\tilde{\textbf{X}}_{\textbf{A}}=\textbf{A}\tilde{\textbf{X}}$, and
$\textbf{z}_{\textbf{A}_{u\mid o}}=\textbf{A}E_{u\mid o}(\textbf{z})$. Then,
we express $\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}$ and
$\textbf{r}_{u\mid o}^{\top}\textbf{M}\textbf{r}_{u\mid o}$ as

$\begin{split}\tilde{\textbf{X}}^{\top}\textbf{M}\tilde{\textbf{X}}&=\tilde{\textbf{X}}^{\top}\textbf{A}^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}\tilde{\textbf{X}}\\\
&=(\textbf{A}\tilde{\textbf{X}})^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}\tilde{\textbf{X}}\\\
&=\tilde{\textbf{X}}_{\textbf{A}}^{\top}\left(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}}(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{\top}\right)^{-1}\tilde{\textbf{X}}_{\textbf{A}}\\\
&=\tilde{\textbf{X}}_{\textbf{A}}^{\top}{(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}}^{\top}(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\tilde{\textbf{X}}_{\textbf{A}}\\\
&=\left((\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\tilde{\textbf{X}}_{\textbf{A}}\right)^{\top}(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\tilde{\textbf{X}}_{\textbf{A}}\\\
&=\textbf{C}^{\top}\textbf{C}.\end{split}$ (S5.10)

and

$\begin{split}\textbf{r}_{u\mid o}^{\top}\textbf{M}\textbf{r}_{u\mid
o}&=(E_{u\mid o}(\textbf{z})-\boldsymbol{\mu})^{\top}\textbf{M}(E_{u\mid
o}(\textbf{z})-\boldsymbol{\mu})\\\ &=E_{u\mid
o}(\textbf{z})^{\top}\textbf{M}E_{u\mid o}(\textbf{z})-2E_{u\mid
o}(\textbf{z})^{\top}\textbf{M}\tilde{\textbf{X}}\boldsymbol{\beta}+\boldsymbol{\beta}^{\top}\textbf{X}^{\top}\textbf{M}\tilde{\textbf{X}}\boldsymbol{\beta}.\end{split}$
(S5.11)

To calculate the quadratic form presented in Equation (S5.11), we need to
calculate $E_{u\mid o}(\textbf{z})^{\top}\textbf{M}E_{u\mid o}(\textbf{z})$,
and $E_{u\mid o}(\textbf{z})^{\top}\textbf{M}\tilde{\textbf{X}}$ as follows:

$\begin{split}E_{u\mid o}(\textbf{z})^{\top}\textbf{M}E_{u\mid
o}(\textbf{z})&=E_{u\mid
o}(\textbf{z})^{\top}\textbf{A}^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}E_{u\mid
o}(\textbf{z})\\\ &=\left(\textbf{A}E_{u\mid
o}(\textbf{z})\right)^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}E_{u\mid
o}(\textbf{z})\\\
&=\left((\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\textbf{z}_{\textbf{A}_{u\mid
o}}\right)^{\top}\left((\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\textbf{z}_{\textbf{A}_{u\mid
o}}\right)\\\ &=\textbf{c}^{\top}\textbf{c}.\end{split}$ (S5.12)

and,

$\begin{split}E_{u\mid
o}(\textbf{z})^{\top}\textbf{M}\tilde{\textbf{X}}&=E_{u\mid
o}(\textbf{z})^{\top}\textbf{A}^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}\tilde{\textbf{X}}\\\
&=(\textbf{A}E_{u\mid
o}(\textbf{z}))^{\top}\left(\textbf{A}\textbf{A}^{\top}+\textbf{I}_{n}\theta\right)^{-1}\textbf{A}\tilde{\textbf{X}}\\\
&=\left((\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\textbf{z}_{\textbf{A}_{u\mid
o}}\right)^{\top}(\textbf{L}_{\textbf{A}\textbf{A}^{\top}+\theta\textbf{I}_{n}})^{-1}\tilde{\textbf{X}}_{\textbf{A}}\\\
&=\textbf{c}^{\top}\textbf{C}.\end{split}$ (S5.13)

##  S6 Computing Information matrix

To obtain variances for the parameters of interest
$\boldsymbol{\beta},\rho,\sigma_{\epsilon}^{2}$ and $\sigma^{2}_{\textbf{y}}$,
the standard procedure requires the calculation of either the expected or the
observed information matrix. They are obtained by first calculating the second
derivatives of the negative of marginal log-likelihood. The negative of
marginal log-likelihood is

$-L_{o}=-\text{log}f(\textbf{z}_{o};\sigma^{2}_{\epsilon},\sigma^{2}_{\textbf{y}},\rho,\boldsymbol{\beta})=\frac{n_{o}}{2}\textrm{log}(2\pi)+\frac{n_{o}}{2}\textrm{log}(\sigma^{2}_{\epsilon})-\frac{1}{2}\textrm{log}|\textbf{V}^{-1}_{oo}|+\frac{1}{2\sigma^{2}_{\epsilon}}\textbf{r}_{o}^{\top}\textbf{V}^{-1}_{oo}\textbf{r}_{o},$
(S6.1)

and for the convenience of notation, we denote the negative marginal log-
likelihood $-L_{o}=\bar{L}_{o}$.

First, we discuss the H-SEM. The second derivative of $\bar{L}_{o}$ with
respect to $\boldsymbol{\beta}$

$\begin{split}\frac{\partial\bar{L}_{o}}{\partial\boldsymbol{\beta}}&=\frac{1}{2\sigma^{2}_{\epsilon}}2(-{\textbf{X}}_{o}^{\top})\textbf{V}_{oo}^{-1}(\textbf{z}_{o}-{\textbf{X}}_{o}\boldsymbol{\beta})\\\
&=-\frac{1}{\sigma_{\epsilon}^{2}}{\textbf{X}}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{r}_{o},\end{split}$
(S6.2)

and since
$E(\textbf{r})=E(\textbf{z}_{o}-{\textbf{X}}_{o}\boldsymbol{\beta})=0$, it is
straightforward to show that

$\displaystyle
E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\rho}\right)=0,$
$\displaystyle\qquad
E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\sigma^{2}_{\epsilon}}\right)=0,$
$\displaystyle\qquad
E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\sigma^{2}_{\textbf{y}}}\right)=0,$

meaning that $\hat{\boldsymbol{\beta}}$ is asymptotically independent from
$\boldsymbol{\zeta}=(\rho,\sigma^{2}_{\epsilon},\sigma^{2}_{\textbf{y}})^{\top}$.
Further,

$E\left(\frac{\partial\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{\top}}\right)=\frac{1}{\sigma_{\epsilon}^{2}}\textbf{X}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{X}_{o}.$
(S6.3)

These results leads to
$\sigma_{\epsilon}^{2}\text{Cov}(\boldsymbol{\hat{\beta}})=\left(\textbf{X}_{o}^{\top}\textbf{V}_{oo}^{-1}\textbf{X}_{o}\right)^{-1}$,
and
$\text{Cov}(\boldsymbol{\hat{\zeta}})=\left[E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\zeta}\partial\boldsymbol{\zeta}^{\top}}\right)\right]^{-1}$.
To calculate $\text{Cov}(\boldsymbol{\hat{\zeta}})$, we recommend utilising
the observed information matrix through numerical differentiation. This
approach helps prevent the need for frequent and expensive inversions. Now, we
consider the H-SAM, which is more complicated because $\textbf{r}_{o}$ depends
on $\rho$. The first derivative of $\bar{L_{o}}$ with respect to
$\boldsymbol{\beta}$ is

$\begin{split}\frac{\partial\bar{L_{o}}}{\partial\boldsymbol{\beta}}&=\frac{1}{2\omega}2(-{(\textbf{A}^{-1}{\textbf{X}})_{o}}^{\top})\textbf{V}_{oo}^{-1}(\textbf{z}_{o}-(\textbf{A}^{-1}{\textbf{X}})_{o}\boldsymbol{\beta})\\\
&=-\frac{1}{\sigma^{2}_{\epsilon}}({\textbf{A}^{-1}{\textbf{X}})_{o}}^{\top}\textbf{V}_{oo}^{-1}\textbf{r}_{o}.\end{split}$
(S6.4)

Then, it is straightforward to show that

$E\left(\frac{\partial\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{\top}}\right)=\frac{1}{\sigma_{\epsilon}^{2}}({\textbf{A}^{-1}{\textbf{X}})_{o}}^{\top}\textbf{V}_{oo}^{-1}({\textbf{A}^{-1}{\textbf{X}})_{o}}.$
(S6.5)

In addition, while
$E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\rho}\right)\neq
0$ , it can be shown that
$E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\sigma^{2}_{\epsilon}}\right)=E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\sigma^{2}_{\textbf{y}}}\right)=0$,
with

$E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\beta}\partial\rho}\right)=\frac{1}{\sigma^{2}_{\epsilon}}\left((\textbf{A}^{-1}\textbf{X})_{o}^{\top}\textbf{V}_{oo}^{-1}\frac{\partial(\textbf{A}^{-1}\textbf{X})_{o}}{\partial\rho}\boldsymbol{\beta}\right),$
(S6.6)

where
$\frac{\partial(\textbf{A}^{-1}\textbf{X})_{o}}{\partial\rho}=(\textbf{A}^{-1}\textbf{W}\textbf{A}^{-1}\textbf{X})_{o}$.
We recommend employing numerical differentiation to compute the corresponding
second derivative for the term
$E\left(\frac{\partial^{2}\bar{L}_{o}}{\partial\boldsymbol{\zeta}\partial\boldsymbol{\zeta}^{\top}}\right)$.
This allows us to derive the observed information matrix without utilising
expensive matrix computations.