[]Yi-JenShih []DavidHarwath

# Interface Design for Self-Supervised Speech Models

###### Abstract

Self-supervised speech (SSL) models have recently become widely adopted for
many downstream speech processing tasks. The general usage pattern is to
employ SSL models as feature extractors, and then train a downstream
prediction head to solve a specific task. However, different layers of SSL
models have been shown to capture different types of information, and the
methods of combining them are not well studied. To this end, we extend the
general framework for SSL model utilization by proposing the interface that
connects the upstream and downstream. Under this view, the dominant technique
of combining features via a layerwise weighted sum can be regarded as a
specific interface. We propose several alternative interface designs and
demonstrate that the weighted sum interface is suboptimal for many tasks. In
particular, we show that a convolutional interface whose depth scales
logarithmically with the depth of the upstream model consistently outperforms
many other interface designs.

###### keywords:

self-supervised speech models, finetuning

##  1 Introduction

Speech processing tasks such as automatic speech recognition have
traditionally relied on supervised learning algorithms that require
transcribed training data. However, because transcriptions are expensive and
time-consuming to collect, self-supervised speech models have recently become
extremely popular. Speech models based on self-supervised learning (SSL) can
be pre-trained with large amounts of unlabeled data, e.g. using a masked
language modeling objective, and subsequently fine-tuned on a relatively small
amount of labeled data for a target downstream task. Recent literature has
proposed various SSL algorithms [1, 2, 3, 4] and many applications have been
built on top of these SSL speech models [5], including Automatic Speech
Recognition [6, 7], Visually Grounded Speech [8, 9, 10, 11], Speech
Segmentation [9, 12] Emotion recognition [13], speaker verification [14] and
so on. Furthermore, there are benchmarks proposed to evaluate speech SSL
models on various downstream tasks [15, 16, 17, 18].

Figure 1: The general framework for utilizing self-supervised speech models
only considers the upstream and downstream model components. We argue that the
interface connecting them should be considered in its own right as a separate
component.($L$ : number of upstream model layers, $T$ : upstream model output
sequence length, $D$ : upstream model feature dimension)

While significant progress has been made in developing algorithms and
applications of speech SSL models, it is still unclear what the best method is
for utilizing these models on downstream tasks. Existing methods for utilizing
speech SSL models fall roughly into 3 categories: 1) Cascade the upstream
model with a task-specific downstream prediction head and fine-tune the entire
model. 2) Freeze the upstream model and take the output of a specific layer as
the input features for a task-specific downstream module, then train the
downstream module. 3) Compute a learnable weighted sum of all layers belonging
to the upstream model and use the result as the input features for a task-
specific downstream module. Typically, only the summation weights and the
downstream module are trained during fine-tuning. Among the three methods,
Fine-tuning the upstream and downstream models together often leads to better
performance, but it is far more computationally expensive and can be unstable
when a small amount of fine-tuning data is available. In contrast, using
features from a single layer of the frozen upstream model is the most
computationally inexpensive, however, it neglects the non-uniform distribution
of different information types encoded across different layers of the upstream
model [19]. The weighted sum can thus be viewed as an intermediate trade-off
between the other two options and is the default configuration used in the
SUPERB [15] benchmark.

However, in this paper, we argue that the weighted sum is still a suboptimal
way of combining information across layers of speech SSL models. Our intuition
is that because the information encoded in each dimension might differ across
each layer of the upstream model, naively summing over the layers dimension-
by-dimension may lead to an information loss in proportion to the degree of
statistical independence between the same dimension across multiple layers.

Built upon this motivation, we hypothesize that there exists a better way of
aggregating the information across different layers and introducing an
Interface module in the SSL framework. The overall framework for utilizing
speech SSL models then becomes a set of three components: the Upstream model,
the Downstream prediction head, and the Interface that bridges between them by
aggregating information across all upstream model layers. Under this
definition, the widely-used weighted sum method is a specific type of
interface, and in this paper, we further propose several alternative interface
designs. We evaluate the performance of each proposed interface across 5
upstream speech SSL models and several downstream tasks in the ML-SUPERB [18]
and SUPERB [15] benchmarks. Among all our proposed interfaces, we show that
Hierarchical Convolution across layers tends to achieve the best overall
performance. We also conduct additional experiments showing that performance
differences among various interface designs are due to their architectures
themselves, and are not merely due to the difference in the number of
trainable parameters. Additionally, we show that the Hierarchical Convolution
interface still outperforms the weighted sum interface when the upstream model
is fine-tuned with the downstream model. The code is publicly available.
111https://github.com/atosystem/SSL_Interface Our contributions in this paper
can be summarized as follows:

  1. 1.

We propose an updated framework for self-supervised speech models that
formally recognizes the interface module

  2. 2.

We introduce several designs for the interface and show that the Hierarchical
Convolution tends to perform the best across many speech processing tasks

  3. 3.

We show that even when the full model is end-to-end fine-tuned, different
interfaces still show performance differences.

##  2 Proposed Interface Methods

Figure 2: Proposed Interface designs. (a): Hierarchical Convolution is applied
over the layer dimension of the upstream model hiddenstates. (b): A learnable
CLS embedding is used for summarizing over layer dimension

###  2.1 Interface Definition

We formalize our proposed overall pipeline for utilizing self-supervised
speech models as a general framework (as shown in Fig. 1): Upstream
$\rightarrow$ Interface $\rightarrow$ Downstream. The Upstream model
$\mathrm{U}\left(\cdot\right):\mathbb{R}^{T^{\prime}}\mapsto\mathbb{R}^{L\times
T\times D}$ stands for a pre-trained, self-supervised speech model which maps
an utterance waveform ${x}\in\mathbb{R}^{T^{\prime}}$ of length $T^{\prime}$
to a hidden representation $\bm{h}$ which consists of $L$ representation
layers each with a feature dimension of $D$. Interface is a function that
aggregates information over the layer dimension $L$ of the upstream output,
$\mathrm{I}\left(\cdot\right):\mathbb{R}^{L\times T\times
D}\mapsto\mathbb{R}^{T\times D}$. Finally, Downstream
$\mathrm{D}\left(\cdot\right)$ takes $\bm{z}=\mathrm{I}\left(\bm{h}\right)$ as
input and outputs predictions for a specific downstream task. Under this
general framework, the weighted sum can be viewed as a specific type of
interface,

$\mathrm{I}_{\mathrm{WS}}\left(\bm{h}\right)=\sum_{l}^{L}{w_{l}\cdot\bm{h}_{l}}$
(1)

where $\bm{w}$ is the per-layer weight vector, which is typically learned
jointly with the downstream model weights.

The weighted sum interface has the benefit of being able to aggregate
information across all $L$ upstream layers with only $L$ learnable parameters
required. On the other hand, we hypothesize that it has potential downsides:
by directly summing layer representations dimension-wise, statistically
independent features from different layers of the upstream model may "collide"
with one another resulting in information loss. As the number of layers in the
upstream model increases, this problem can become more severe. To further
investigate this, we propose several alternative interface designs that
attempt to avoid information collision across multiple upstream layers while
remaining relatively lightweight and also keeping the output dimension roughly
the same as the upstream model's dimension.

###  2.2 Proposed Interfaces

####  2.2.1 Grouped Weighted Sums

An intuitive way to reduce the information collision of a single weighted sum
is instead to use multiple weighted sums, each of which only has access to a
subset of the upstream layers. Specifically, after conducting a weighted sum
within each subset of upstream layers, we concat the results of each group
along the feature dimension. Then we use a learnable projection to project the
feature dimension to that of the downstream model.

####  2.2.2 Concatenation + Learnable Projection

Another type of interface concatenates all layer representations $\bm{h}$
along the feature dimension and uses a learnable projection layer to
automatically decide the features that are beneficial to the downstream model.
Specifically, we reshape $\bm{h}$ from $\left(L,T,D\right)$ to
$\left(T,L\times D\right)$, then project $\left(T,L\times D\right)$ back to
the original upstream model dimension $\left(T,D\right)$.

####  2.2.3 Hierarchical Convolution over Layers

Even though the feature manifold differs among each layer of the upstream
model, neighboring layers can still possess similar feature distributions due
to the residual connection in the transformer layers [20]. Motivated by the
local structure of the hidden representation across the layer dimension, we
apply 1D convolutions over the layer dimension to aggregate the information
across layers. In practice, an interface design should be adaptable to
upstream models with an arbitrary number of layers. Hence, in our paper, we
fix the kernel size=$5$, stride=$3$ for our convolution layer and we stack
$\lfloor\mathrm{log}_{3}{L}\rfloor$ identical convolutional layers, which has
the effect of collapsing the output of all models layers down to a single
vector at the output of the convolutional stack (as shown in Fig. 2 (a)). We
leave further convolutional hyperparameter optimization for future study.

####  2.2.4 CLS Pooling over layer dimension

We can use an attention module to aggregate information across layers rather
than using convolutions. Specifically, at each timestep we view all the layers
of the upstream model as a sequence of tokens and concatenate a learnable CLS
vector [21] to this sequence. A Transformer layer is then used to aggregate
information across all layers at the same time step into the CLS token (as
shown in Fig. 2 (b)). A key difference compared to our other proposed
interface designs is that CLS Pooling is data-dependent, while other methods
are only task-dependent.

####  2.2.5 PCA + Concatenation

Besides the above parametric methods, we also investigated layer-wise
dimensionality reduction via Principal Component Analysis (PCA), followed by
concatenation across layers. Using the HuBERT Base model as an example
upstream model, we take only the top 60 principal components of each
Transformer layer and concatenate them together, resulting in $60\times
13=780\approx 768$ dimensions, which is roughly the same as the upstream model
feature dimension. Practically, the PCA transform is learned from each
downstream task dataset respectively.

##  3 Experiments

Table 1: Comparison between different proposed interfaces with HuBERT Base as the upstream model, evaluated on the ML-SUPERB Monolingual track. ``GroupWS'' indicates Grouped Weighted Sums and ``Hierarchical Conv.'' stands for Hierarchical Convolution over layers. The metric is WER/CER. Interface | Params | Mono-10min | Mono-1hr
---|---|---|---
Weighted Sum | 13 | 42.85 | 35.15
GroupWS(#Groups=2) | 1.2M | 41.84 | 34.47
GroupWS(#Groups=3) | 1.8M | 43.08 | 33.96
GroupWS(#Groups=4) | 2.4M | 42.52 | 33.99
Concat $+$ Proj | 7.7M | 43.20 | 34.26
PCA $+$ Concat | 0 | 45.16 | 36.58
Hierarchical Conv. | 4.4M | 41.51 | 33.88
CLS Pooling | 5.5M | 48.36 | 33.92
Table 2: Performance comparison of proposed interfaces using different upstream models, evaluated on both ML-SUPERB and SUPERB. The metric for each task: ``Mono-1hr'' and ``Multi-1hr'': WER/CER, ``LID-1hr'',``ER'',``IC'', ``SV'': accuracy, ``PR'': PER.  Upstream | Interface | ML-SUPERB | SUPERB
---|---|---|---
Mono-1hr ($\downarrow$) | Multi-1hr ($\downarrow$) | LID-1hr ($\uparrow$) | ER ($\uparrow$) | IC ($\uparrow$) | SV ($\downarrow$) | PR ($\downarrow$)
HuBERT Base | Weighted Sum | 35.1 | 24.4 | 84.6 | 64.92 | 98.34 | 5.11 | 5.41
Hierarchical Conv. | 33.9 | 24.0 | 81.7 | 68.49 | 99.45 | 5.62 | 3.07
CLS Pooling | 33.9 | 24.3 | 1.2 | 62.73 | 99.45 | 49.91 | 2.93
HuBERT Large | Weighted Sum | 32.3 | 22.3 | 64.4 | 67.62 | 98.76 | 5.98 | 3.53
Hierarchical Conv. | 30.0 | 21.4 | 89.2 | 72.44 | 99.53 | 6.03 | 1.76
CLS Pooling | 31.1 | 21.5 | 87.0 | 68.72 | 99.53 | 6.34 | 1.64
WavLM Base | Weighted Sum | 34.2 | 24.3 | 84.1 | 65.94 | 98.63 | 4.69 | 4.84
Hierarchical Conv. | 32.4 | 23.6 | 67.9 | 68.57 | 99.53 | 5.48 | 3.06
CLS Pooling | 33.8 | 23.8 | 13.2 | 63.60 | 99.50 | 6.03 | 2.92
WavLM Large | Weighted Sum | 30.1 | 20.8 | 71.7 | 70.62 | 99.31 | 3.77 | 3.06
Hierarchical Conv. | 28.0 | 19.4 | 90.6 | 74.95 | 99.71 | 5.20 | 1.72
CLS Pooling | 29.7 | 18.8 | 88.5 | 70.45 | 99.63 | 4.27 | 1.73
XLSR-53 | Weighted Sum | 35.1 | 20.2 | 76.6 | 66.34 | 95.62 | 6.45 | 4.50
Hierarchical Conv. | 30.6 | 19.8 | 80.4 | 72.01 | 99.55 | 5.63 | 2.69
CLS Pooling | 31.0 | 80.7 | 4.5 | 65.10 | 99.47 | 12.95 | 5.27

###  3.1 Experimental Setup

To get a holistic performance comparison of different interface designs, we
choose 5 SSL models and 2 speech SSL benchmarks for evaluation. Specifically,
we pick HubERT Base and Large [1], WavLM Base and Large [2], and XLSR-53 [3]
for our upstream models. For downstream tasks, we test all proposed interfaces
on ML-SUPERB [18] because it is currently less saturated than many SUPERB
tasks. ML-SUPERB contains both 10min and 1hr training sets for each language.
In the ML-SUPERB Monolingual track, the overall reported score is the average
CER across 13 languages (the model for each language is trained individually).
In the multilingual track, there are two tasks: ASR and Language
Identification (LID). In these tasks, a model must be trained on 143 languages
simultaneously.

In addition to multilingual ASR tasks, we also test the generalization of our
interface designs on other common speech processing tasks. Hence, we select 4
tasks from SUPERB: Intent Classification (IC), Phoneme Recognition (PR),
Emotion Recognition (PR), and Speaker Verification (SV).

To reduce the computational cost of running a large number of potential
experiments, we first test all proposed interfaces on ML-SUPERB Monolingual
Track with HuBERT Base. Then we pick the top 2 most promising interface
designs and run them on all 5 upstream models and both ML-SUPERB and SUPERB.
For both SUPERB and ML-SUPERB, we follow the default training configurations
as far as possible. However, for the Multilingual track on ML-SUPERB, we find
that the default downstream model is not optimal, and performance can be
easily increased by simply adding more layers. In order to disentangle the
effect of a larger downstream model from the interface itself, we do a
hyperparameter search on the best number of layers in the downstream model so
that it can achieve the lowest CER in terms of the same frozen upstream model.
From our empirical experiment, we found that 65 is the optimal layer number
for Multilingual ASR (about 56M params).

###  3.2 Interface Comparison on ML-SUPERB Monolingual ASR with HuBERT Base

The results for our pilot experiment using HuBERT Base on the ML-SUPERB
Monolingual task are shown in Table 1. Firstly, we notice that there is a
significant performance difference between all proposed interfaces relative to
the baseline weighted sum method. While several of the proposed interfaces
improve over the baseline weighted sum, we see that the Hierarchical
Convolution achieves the overall best interface performance on both the 10min
and 1hr datasets. The CLS Pooling interface ranks second on the 1hr task, but
performs poorly on 10min, indicating that this method likely requires more
finetuning data to reach its potential. In light of their promising results on
the 1hr ML-SUPERB monolingual task, we choose Hierarchical Convolution and CLS
Pooling for further experimentation.

###  3.3 Interface Comparison on ML-SUPERB and SUPERB

The results for the full suite of tasks across ML-SUPERB and SUPERB using all
upstream models are shown in Table 2. Again, we observe that the Hierarchical
Convolution interface tends to consistently outperform the baseline weighted
sum as well as the CLS Pooling interface. The performance differences are
often substantial: for example, simply replacing the weighted sum interface
with the Hierarchical Convolution interface reduces HuBERT Base's PER on the
SUPERB phone recognition task from 5.41 to 2.93, which is lower than the 3.53
PER achieved by HuBERT Large with the weighted sum interface. We even find
that the CLS Pooling interface with the HuBERT Large model establishes a new
SotA on the SUPERB leaderboard for phone recognition at 1.64 PER, beating the
WavLM model's 3.06 PER when using the weighted sum.

Generally, the Hierarchical Convolution interface outperform the baseline
weighted sum on all tasks except for Speaker Verification with HuBERT and
WavLM models, although XLSR-53 with Hierarchical Convolution outperforms the
weighted sum. This might imply that speaker-related information is encoded
similarly across different layers in HuBERT and WavLM, causing our proposed
interfaces to overfit.

Interestingly, in both benchmarks, we also observe a larger improvement over
the basic weighted sum in Large models compared to Base models for HuBERT and
WavLM when Hierarchical Convolution is applied, which supports our hypothesis
that the basic weighted sum averages out more information as the number of
layers increases. In particular, on the LID task, we see that with Base
upstream models the Hierarchical Convolution interface underperforms the basic
weighted sum, but when scaling to the Large models (HuBERT Large, WavLM Large,
and XLSR-53) the Hierarchical Convolution instead outperforms the basic
weighted sum.

###  3.4 Interface vs. Larger Downstream Models

Table 3: Ablation experiments on trainable parameters. ``WS w/ Large DS'' indicates Weighted Sum with Large Downstream model. Upstream | Interface | ML-SUPERB (WER/CER)
---|---|---
Mono-1hr ($\downarrow$) |  Multi-1hr ($\downarrow$)
HuBERT Base | Hierarchical Conv. | 33.9 | 24.0
WS w/ Large DS | 35.2 | 24.8
HuBERT Large | Hierarchical Conv. | 30.0 | 21.4
WS w/ Large DS | 32.9 | 22.1
WavLM Base | Hierarchical Conv. | 32.4 | 23.6
WS w/ Large DS | 34.4 | 24.3
WavLM Large | Hierarchical Conv. | 28.0 | 19.4
WS w/ Large DS | 30.5 | 20.8
XLSR-53 | Hierarchical Conv. | 30.6 | 19.8
WS w/ Large DS | 35.4 | 20.0

Although there are substantial improvements using the Hierarchical Convolution
interface, it is not immediately clear whether the performance gain comes from
the interface design or simply the additional trainable parameters. To this
end, we increase the size of the downstream model for Monolingual and
Multilingual ASR on ML-SUPERB by an amount roughly the same as the number of
parameters in the Hierarchical Convolution interface (denoted as ``WS+Large
DS'').

As shown in Table 3, under all settings, weighted sum with a large downstream
model does not lead to better performance than the Hierarchical Convolution
interface combined with a smaller downstream model. This implies the role of
the interface is not simply to add more learnable parameters to the downstream
model, but to actually better aggregate information across the layers of the
upstream model. Our findings not only resonate with [22] that larger
downstream models increase the performance but we also show that a better
interface design can get further improvements.

###  3.5 Performance under Fine-tuning setting

We motivated the introduction of the interface module under the assumption
that the upstream model would typically remain frozen during the entire
downstream task fine-tuning. However, we wish to empirically verify whether
different interfaces can still offer improvements with end-to-end fine-tuning.
Practically, we follow the same pipeline as in our previous experiments except
the upstream model is trainable. We evaluate on ML-SUPERB Monolingual track 1h
using HuBERT Base.

Table 4: Comparison of Weighted Sum and Hierarchical under end-to-end fine-tuning. The metric is WER/CER. Interface |  ML-SUPERB Mono-1hr ($\downarrow$)
---|---
Freeze Upstream | Fine-tuned
Weighted Sum | 35.1 | 31.5
Hierarchical Conv. | 33.9 | 31.1
WS w/ Large DS | 35.2 | 31.6

As shown in Table 4, end-to-end fine-tuning shortens the gap between different
interface designs. However, the Hierarchical Convolution interface still
outperforms the weighted sum baseline. These results indicate that the
interface design is beneficial regardless of frozen or trainable upstream
model.

##  4 Conclusion

In this work, we introduced the interface module into the conceptual framework
for utilizing self-supervised speech models. We also proposed several new
interface designs. From our experiments on multiple benchmarks and upstream
models, we showed that the Hierarchical Convolution interface tends to
outperform our other proposed designs, as well as the weighted sum method that
is currently standard practice in the field. Finally, our ablation experiments
substantiate that the role of the interface in optimally aggregating
information across layers of the upstream model is more important than the
amount of trainable parameters during funetuning.

##  5 Acknowledgements

This work is supported in part by the National Science Foundation under Grant
No. 2238605.

## References

  * [1] W.-N. Hsu, B. Bolte, Y.-H. H. Tsai, K. Lakhotia, R. Salakhutdinov, and A. Mohamed, ``Hubert: Self-supervised speech representation learning by masked prediction of hidden units,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 29, pp. 3451–3460, 2021.
  * [2] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao, J. Wu, L. Zhou, S. Ren, Y. Qian, Y. Qian, M. Zeng, and F. Wei, ``Wavlm: Large-scale self-supervised pre-training for full stack speech processing,'' _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, pp. 1505–1518, 2021.
  * [3] A. Conneau, K. Khandelwal, N. Goyal, V. Chaudhary, G. Wenzek, F. Guzmán, E. Grave, M. Ott, L. Zettlemoyer, and V. Stoyanov, ``Unsupervised cross-lingual representation learning at scale,'' in _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_. Association for Computational Linguistics, Jul. 2020, pp. 8440–8451.
  * [4] S. Schneider, A. Baevski, R. Collobert, and M. Auli, ``wav2vec: Unsupervised pre-training for speech recognition,'' in _Interspeech_ , 2019.
  * [5] A. Mohamed, H.-y. Lee, L. Borgholt, J. D. Havtorn, J. Edin, C. Igel, K. Kirchhoff, S.-W. Li, K. Livescu, L. Maaløe _et al._ , ``Self-supervised speech representation learning: A review,'' _IEEE Journal of Selected Topics in Signal Processing_ , 2022.
  * [6] A. Baevski, W.-N. Hsu, A. Conneau, and M. Auli, ``Unsupervised speech recognition,'' _NeurIPS_ , 2021.
  * [7] A. H. Liu, W.-N. Hsu, M. Auli, and A. Baevski, ``Towards end-to-end unsupervised speech recognition,'' in _2022 IEEE Spoken Language Technology Workshop (SLT)_ , 2023, pp. 221–228.
  * [8] P. Peng and D. Harwath, ``Fast-slow transformer for visually grounding speech,'' in _ICASSP_ , 2022.
  * [9] ——, ``Word discovery in visually grounded, self-supervised speech models,'' in _Interspeech_ , 2022.
  * [10] Y.-J. Shih, H.-F. Wang, H.-J. Chang, L. Berry, H.-y. Lee, and D. Harwath, ``Speechclip: Integrating speech with pre-trained vision and language model,'' in _2022 IEEE Spoken Language Technology Workshop (SLT)_ , 2023, pp. 715–722.
  * [11] L. Berry, Y.-J. Shih, H.-F. Wang, H.-J. Chang, H.-Y. Lee, and D. Harwath, ``M-speechclip: Leveraging large-scale, pre-trained models for multilingual speech to image retrieval,'' in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5.
  * [12] L. Strgar and D. Harwath, ``Phoneme segmentation using self-supervised speech models,'' in _2022 IEEE Spoken Language Technology Workshop (SLT)_ , 2023, pp. 1067–1073.
  * [13] E. Morais, R. Hoory, W. Zhu, I. Gat, M. Damasceno, and H. Aronowitz, ``Speech emotion recognition using self-supervised features,'' in _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2022, pp. 6922–6926.
  * [14] J. Peng, T. Stafylakis, R. Gu, O. Plchot, L. Mošner, L. Burget, and J. Černocký, ``Parameter-efficient transfer learning of pre-trained transformer models for speaker verification using adapters,'' in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5.
  * [15] S. Yang, P.-H. Chi, Y.-S. Chuang, C.-I. J. Lai, K. Lakhotia, Y. Y. Lin, A. T. Liu, J. Shi, X. Chang, G.-T. Lin, T.-H. Huang, W.-C. Tseng, K. tik Lee, D.-R. Liu, Z. Huang, S. Dong, S.-W. Li, S. Watanabe, A. Mohamed, and H. yi Lee, ``SUPERB: Speech Processing Universal PERformance Benchmark,'' in _Proc. Interspeech 2021_ , 2021, pp. 1194–1198.
  * [16] H.-S. Tsai, H.-J. Chang, W.-C. Huang, Z. Huang, K. Lakhotia, S.-w. Yang, S. Dong, A. Liu, C.-I. Lai, J. Shi, X. Chang, P. Hall, H.-J. Chen, S.-W. Li, S. Watanabe, A. Mohamed, and H.-y. Lee, ``SUPERB-SG: Enhanced speech processing universal PERformance benchmark for semantic and generative capabilities,'' in _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , S. Muresan, P. Nakov, and A. Villavicencio, Eds. Dublin, Ireland: Association for Computational Linguistics, May 2022, pp. 8479–8492. [Online]. Available: https://aclanthology.org/2022.acl-long.580
  * [17] S. Shon, S. Arora, C.-J. Lin, A. Pasad, F. Wu, R. S. Sharma, W.-L. Wu, H.-y. Lee, K. Livescu, and S. Watanabe, ``SLUE phase-2: A benchmark suite of diverse spoken language understanding tasks,'' in _Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_. Association for Computational Linguistics, Jul. 2023, pp. 8906–8937. [Online]. Available: https://aclanthology.org/2023.acl-long.496
  * [18] J. Shi, D. Berrebbi, W. Chen, E.-P. Hu, W.-P. Huang, H.-L. Chung, X. Chang, S.-W. Li, A. Mohamed, H. yi Lee, and S. Watanabe, ``ML-SUPERB: Multilingual Speech Universal PERformance Benchmark,'' in _Proc. INTERSPEECH 2023_ , 2023, pp. 884–888.
  * [19] A. Pasad, B. Shi, and K. Livescu, ``Comparative layer-wise analysis of self-supervised speech models,'' in _ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023, pp. 1–5.
  * [20] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, ``Attention is all you need,'' in _NeurIPS_ , 2017.
  * [21] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, ``BERT: Pre-training of deep bidirectional transformers for language understanding,'' in _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , J. Burstein, C. Doran, and T. Solorio, Eds. Minneapolis, Minnesota: Association for Computational Linguistics, Jun. 2019, pp. 4171–4186. [Online]. Available: https://aclanthology.org/N19-1423
  * [22] S. Zaiem, Y. Kemiche, T. Parcollet, S. Essid, and M. Ravanelli, ``Speech Self-Supervised Representation Benchmarking: Are We Doing it Right?'' in _Proc. INTERSPEECH 2023_ , 2023, pp. 2873–2877.


# Efficient Estimation For The Joint Model of Survival and Longitudinal Data

###### Abstract

In survival studies it is important to record the values of key longitudinal
covariates until the occurrence of event of a subject. For this reason, it is
essential to study the association between longitudinal and time-to-event
outcomes using the joint model. In this paper, profile likelihood approach has
been used to estimate the cumulative hazard and regression parameters from the
joint model of survival and longitudinal data. Moreover, we show the
asymptotic normality of the profile likelihood estimator via asymptotic
expansion of the profile likelihood and obtain the explicit form of the
variance estimator. The estimators of the regression parameters from the joint
model are shown to be semi-parametric efficient. The numerical result of our
proposed method is shown by the ‘aids’ data (from JM R-package) where we have
computed the standard errors of the estimated parameters from the efficient
score function and efficient information matrix.

Cox PH model,

Efficient score Function,

EM algorithm,

Implicitly function,

Joint model,

Longitudinal data,

Profile likelihood,

Survival data,

###### keywords:

Khandoker Akib Mohammad∗, Yuichi Hirose, Budhi Surya and Yuan Yao

Victoria University of Wellington

##  1 Introduction

In survival analysis, joint model of repeated measurements and time-to-event
data has attracted attention in the last decade. In survival analysis, many
clinical studies often deal with two types of outcome such as longitudinal
response measurements and the event time of interest (e.g. death, development
of disease etc.) [Rizopoulos, 2010]. In classical survival analysis, the
longitudinal outcome and the event outcome can be modelled separately by using
mixed effect model and Cox PH model respectively. However, to see the effect
of longitudinal outcome as a time-dependent covariate on event outcome,
Faucett and Thomas (1996) and Wulfsohn and Tsiatis (1997) proposed the joint
models of longitudinal and survival outcomes where the longitudinal outcome
and the event outcome were modelled jointly rather than separately. In
clinical settings, the joint model of longitudinal and survival data has been
widely used for various types of medical studies (Brown and Ibrahim, 2003;
Ratcliffe, Guo and Have, 2004; Heish, Tseng and Wang, 2006; Rizopoulos, 2010;
Hanson, Branscum and Johnson, 2011; Viviani, Alfo and Rizopoulos, 2014)

The efficiency and asymptotic distribution of semiparametric maximum
likelihood estimator have been studied for the joint model by Zeng and Cai
(2005). Moreover, Hanson et al. (2011), Chi and Ibrahim (2006), Brown and
Ibrahim (2003), Xu and Zeger (2001), and Wang and Taylor (2001) considered the
Bayesian estimation approach of joint models using MCMC technique. Tsiatis and
Davidian (2001) developed a set of unbiased estimating equations that yields
consistent and asymptotically normal estimators where random effects have
treated as nuisance parameters in their proposed conditional score approach.
Though the non-parametric maximum likelihood estimation has been proposed by
Wulfsohn and Tsiatis (1997) but they did not address the computation with
implicit function in the profile likelihood estimation. Rizopoulos (2010)
developed a R package (JM) to fit the joint model of longitudinal and survival
data wher Gaussian quadrature rules have been used for the numerical
integration with respect to the random effects. However, computational
complexity increases with the dimension of the random-effects. For this
reason, Rizopoulos (2012) proposed a pseudo-adaptive Gauss–Hermite quadrature
rule which is considerably faster than the standard Gauss–Hermite rule. In JM
package, Rizopoulos (2010) used the aids data where the standard errors of the
estimated parameters have been calculated by using bootstrap methods.

In this paper, we profile out the baseline hazard function from the joint
model and plugged the estimator in the likelihood function. However the
problem is that the estimator of the baseline hazard function is an implicit
function (Rizopoulos, 2012). We can solve the difficulty without
differentiating the implicit function to show the asymptotic normality of the
estimator (Theorem-3 and Theorem-4 in Section 3.3). This approach is
alternative to the methodologies where the asymptotic normality of profile
likelihood estimator has been studied (Hirose, 2011, 2016; Murphy and Vaart,
2016). By using the asymptotic expansion of the likelihood, we have found the
explicit form of the efficient score function and established the asymptotic
normality of the profile likelihood estimator. Hence we got the explicit form
of the estimate of variance for the profile likelihood estimator. These
results can be used in computation to calculate the variance of the profile
likelihood estimator which is illustrated in the numerical example
(Section-4). For the actual computation, we have used the same data (aids from
JM R-package) and computed the standard errors of the estimated parameters
from the efficient score function and efficient information matrix.

This paper is organized as follows. A brief discussion on joint model of
longitudinal and survival data has been given in Section-2. In Section-3, we
describe the estimation procedure and theorems which are used to show that the
profile likelihood estimators are consistent and asymptotically normal.
Results obtained from the profile expansion of the joint model are shown in
Section-4. This paper concludes in Section-5 with a short discussion.

##  2 Joint Model of Survival and Longitudinal Data

For each subject $i=1,\ldots,n$, we observe response $y_{i}(t)$, fixed effects
covariate vector $x_{i}(t)$ and random effects covariate vector $z_{i}(t)$.
Now we assume a linear mixed model for the process as

$\displaystyle y_{i}(t)$ $\displaystyle=$ $\displaystyle
m_{i}(t)+\varepsilon_{i}$

where $m_{i}(t)=x^{\prime}_{i}(t)\beta+z^{\prime}_{i}(t)b_{i}$ is the
longidudinal process with the fixed effects $\beta$ and the random effects
$b_{i}$. We assume the error $\varepsilon_{i}\sim N(0,\sigma^{2}I_{n_{i}})$
and the random effects $b_{i}\sim N(0,D)$ are independent, where $I_{n_{i}}$
denotes the $n_{i}$-dimensionality identity matrix and $D$ is the covariance
matrix.

In addition, for the survival data, let us assume $T^{*}_{i}$ and $C_{i}$ be
the survival and censoring times, respectively. Moreover, $T_{i}$ be the
observed event time which can be defined as $T_{i}=\min(T^{*}_{i},C_{i})$ and
$\delta_{i}=I(T^{*}_{i}\leq C_{i})$ be the censoring indicator. We observe
$(T_{i},\delta_{i},w_{i})$, where $w_{i}$ be a vector of baseline covariates.
We assume there is no measurements after the event or the censorship:
$t_{in_{i}}\leq T_{i}$.

To measure the strength of association between $m_{i}(t)$ and the risk of the
event, we can express the relative risk models as

$\displaystyle\lambda_{i}(t|b_{i})=\lambda_{0}(t)\exp\left\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(t)\right\\}.$
(2.1)

where $\lambda_{0}(t)$ is the baseline hazard function and $\gamma$ is the
vector of regression coefficients corresponding to $w_{i}$. Similarly $\alpha$
represents the effect of the underlying longitudinal outcome to the risk of
the event.

Let $\Lambda(t)=\int_{0}^{t}\lambda_{0}(s)ds$ be the cumulative baseline
hazard function. Then Survival function can be written as

$S_{i}(t|b_{i})=\exp\left(-\int_{0}^{t}\lambda_{0}(s)\exp\left\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(s)\right\\}ds\right),$

##  3 Estimation

To construct the joint model, we assume that random effects $b_{i}$ underlies
both the longitudinal and event outcomes (conditional independence). The joint
likelihood can be written as

$L(\theta,\Lambda)=\prod_{i=1}^{n}\int_{b_{i}}p(T_{i},\delta_{i},y_{i},b_{i};\theta,\Lambda)db_{i}=\prod_{i=1}^{n}\int_{b_{i}}p(T_{i},\delta_{i}|b_{i};\theta,\Lambda)p(y_{i}|b_{i};\theta)p(b_{i};\theta)db_{i}$

where

$\displaystyle p(T_{i},\delta_{i}|b_{i};\theta,\Lambda)$ $\displaystyle=$
$\displaystyle\left[\lambda_{0}(T_{i})\exp\left\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(T_{i})\right\\}\right]^{\delta_{i}}$
(3.1)
$\displaystyle\times\exp\left(-\int_{0}^{T_{i}}\exp\left\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(s)\right\\}d\Lambda(s)\right)$
$\displaystyle p(y_{i}|b_{i};\theta)$ $\displaystyle=$
$\displaystyle(2\pi\sigma^{2})^{-n_{i}/2}\exp\left\\{-\frac{\sum_{j=1}^{n_{i}}[y_{i}(t_{ij})-m_{i}(t_{ij})]^{2}}{2\sigma^{2}}\right\\},$
(3.2)

and

$\displaystyle p(b_{i}|\theta)$ $\displaystyle=$
$\displaystyle(2\pi)^{-r/2}\det(D)^{-1/2}\exp\left\\{-\frac{b_{i}^{\prime}D^{-1}b_{i}}{2}\right\\},$
(3.3)

where $r$ is the dimension of $b_{i}$ and
$\theta=(\alpha^{\prime},\beta^{\prime},\gamma^{\prime},\sigma^{\prime},\textrm{vec}(D)^{\prime})^{\prime}$
be the all parameters in the joint model except the cumulative baseline hazard
function $\Lambda$. Moreover, “$\textrm{vec}(D)$” is the vectorization of the
matix $D$. In our joint model, $\theta$ is the main parameter of interest
whereas $\Lambda$ is the nuisance parameter.

###  3.1 EM Algorithm

Let us define the complete data by
$(t_{i},\delta_{i},y_{i},b_{i}),~{}i=1,...,n$ which includes the observed data
and unobserved $b_{i}$ (treating $b_{i}$ as missing data). The choice for
using EM algorithm is justified by the fact that the model depends on a latent
variable, $b_{i}$. Moreover, the aim of EM algorithm is to maximize observed
data likelihood from a complete data likelihood. So the complete data log-
likelihood function can be written as

$\log
L_{c}(\theta,\Lambda;b)=\sum_{i=1}^{n}\log\bigg{\\{}p(T_{i},\delta_{i}|b_{i};\theta,\Lambda)p(y_{i}|b_{i};\theta)p(b_{i};\theta)\bigg{\\}}.$

Now the expected complete data log-likelihood under $p(b|T,\delta,y)$ will be

$\displaystyle E_{b}[\log L_{c}(\theta,\Lambda;b)]={}$
$\displaystyle\sum_{i=1}^{n}\int_{b_{i}}\log\bigg{\\{}p(T_{i},\delta_{i}|b_{i};\theta,\Lambda)p(y_{i}|b_{i};\theta)p(b_{i};\theta)\bigg{\\}}p(b_{i}|T_{i},\delta_{i},y_{i};\theta,\Lambda)db_{i}$
(3.4) $\displaystyle={}$
$\displaystyle\sum_{i=1}^{n}\bigg{\\{}E_{b}\big{[}\log
p(T_{i},\delta_{i}|b_{i};\theta,\Lambda)\big{]}+E_{b}[\log
p(y_{i}|b_{i};\theta)]+E_{b}[\log p(b_{i};\theta)]\bigg{\\}}$

where $E_{b}$ is the expectation with respect to $b$ under $p(b|T,\delta,y)$
which can be expressed as

$p(b|T,\delta,y;\theta,\Lambda)=\frac{p(T,\delta|b;\theta,\Lambda)p(y|b;\theta)p(b;\theta)}{\int_{b}p(T,\delta|b;\theta,\Lambda)p(y|b;\theta)p(b;\theta)db}\propto
p(T,\delta|b;\theta,\Lambda)p(y|b;\theta)p(b;\theta)$

### Baseline Hazard Estimation

Before starting the EM algorithm, we profile out the baseline hazard function
$\lambda_{0}(t)$ using NPMLE (non-parametric maximum likelihood estimation).
Without loss of generality we assume observed times $t_{i}$ are ordered:
$t_{1}<t_{2}<\ldots<t_{n}$. We assume $\lambda_{i}$ be the hazard at time
$T_{i}$ and treat the survival time as discrete. The cumulative bseline hazard
function is then $\Lambda_{0}(t_{i})=\sum_{j=1}^{n}\lambda_{j}1\\{t_{j}\leq
t_{i}\\}$.

Now the survival part of equation (3.4) can be written as

$\sum_{i=1}^{n}E_{b}\bigg{[}\delta_{i}\big{\\{}\log\lambda_{i}+\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(t)\big{\\}}-\exp\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(t)\\}\sum_{j=1}^{n}\lambda_{j}1\\{t_{j}\leq
t_{i}\\}\bigg{]}.$ (3.5)

Now from the derivative with respect to $\lambda_{k}$, we get

$\hat{\lambda}_{k}(t)=\frac{\delta_{k}}{\sum_{l=1}^{n}1\\{t_{k}\leq
t_{l}\\}E_{b}[\exp\\{\gamma^{\prime}w_{l}+\alpha^{\prime}m_{l}(t_{k})\\}]}.$

So the estimate of the baseline cumulative hazard will be

$\hat{\Lambda}(t)=\sum_{i=1}^{n}\frac{\delta_{i}1\\{t_{i}\leq
t\\}}{\sum_{l=1}^{n}1\\{t\leq
t_{l}\\}E_{b}[\exp\\{\gamma^{\prime}w_{l}+\alpha^{\prime}m_{l}(t)\\}]}.$ (3.6)

### The E-step

In the E-step, we use the current parameter estimates $\theta$ to find the
expected values of $b$ under $p(b|T,\delta,y;\theta,\hat{\Lambda})$:

$p(b|T,\delta,y;\theta,\hat{\Lambda})\propto
p(T,\delta|b;\theta,\hat{\Lambda})p(y|b;\theta)p(b;\theta)$ (3.7)

### The M-step

By replacing $\lambda$ with $\hat{\lambda}(\theta)$, we maximize the equation
(3.4)

$\displaystyle E_{b}[\log L_{c}(\theta,\hat{\Lambda};b)]={}$
$\displaystyle\sum_{i=1}^{n}\bigg{\\{}E_{b}\big{[}\log
p(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda})\big{]}+E_{b}[\log
p(y_{i}|b_{i};\theta)]+E_{b}[\log p(b_{i};\theta)]\bigg{\\}}$ (3.8)

with respect to $\theta$ to obtain $\hat{\theta}$ respectively. The estimated
parameters from the M-step are returned into E-step until the values of
$\hat{\theta}$ converge.

###  3.2 Score Functions

An estimator of the baseline cumulative hazard function in the counting
process notation (fleming2011counting) can be written from equation (3.6) as

$\hat{\Lambda}(t)=\int_{0}^{t}\frac{\sum_{i=1}^{n}dN_{i}(u)}{\sum_{i=1}^{n}Y_{i}(u)E_{b}[\exp\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(u)\\}]},$
(3.9)

where $N_{i}(u)=1\\{T_{i}\leq u,\delta=1\\}$ and $Y_{i}(u)=1\\{T_{i}\geq
u\\}.$

Let us denote $E_{F_{n}}f=\int fdF_{n}$. Then $\hat{\Lambda}(t)$ can be
expressed as

$\hat{\Lambda}_{\theta,F_{n}}(t)=\int_{0}^{t}\frac{E_{F_{n}}dN(u)}{E_{F_{n}}Y(u)E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\\}]}.$
(3.10)

Now from (3.8), the log-profile likelihood can be written as

$\displaystyle E_{b}[\log L_{c}(\theta,\hat{\Lambda}_{\theta,F_{n}};b)]={}$
$\displaystyle\sum_{i=1}^{n}\bigg{\\{}E_{b}\big{[}\log
p(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F_{n}})\big{]}+E_{b}[\log
p(y_{i}|b_{i};\theta)]+E_{b}[\log p(b_{i};\theta)]\bigg{\\}}$ (3.11)

where $\log p(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F_{n}})$,
$\log p(y_{i}|b_{i};\theta)$ and $\log p(b_{i};\theta)$ are the log-profile
likelihood functions (for one observation) for Cox PH, longitudinal and random
effect components respectively. Now we can express the components as

$\displaystyle E_{b}\big{[}\log
p(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F_{n}})\big{]}={}$
$\displaystyle
E_{b}\bigg{[}\delta_{i}\bigg{\\{}\log\frac{E_{F_{n}}dN(T_{i})}{E_{F_{n}}Y(T_{i})E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(T_{i})\\}]}+\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(T_{i})\bigg{\\}}$
(3.12)
$\displaystyle-\exp\bigg{(}\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(T_{i})\bigg{)}\int_{0}^{T_{i}}\frac{E_{F_{n}}dN(u)}{E_{F_{n}}Y(u)E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\\}]}\bigg{]},$

$E_{b}[\log
p(y_{i}|b_{i};\theta)]=E_{b}\bigg{[}-\frac{1}{2}\bigg{\\{}n_{i}\log(2\pi)+\log(\sigma^{2})+(y_{i}-x_{i}\beta-
z_{i}b_{i})^{\prime}(\sigma^{2})^{-1}(y_{i}-x_{i}\beta-
z_{i}b_{i})^{\prime}\bigg{\\}}\bigg{]},$ (3.13)

and

$E_{b}[\log
p(b_{i};\theta)]=E_{b}\bigg{[}-\frac{1}{2}\bigg{\\{}r\log(2\pi)+\log|D|+b_{i}^{\prime}D^{-1}b_{i}\bigg{\\}}\bigg{]}.$
(3.14)

The score functions for the profile likelihood are

$S(T_{i},\delta_{i},y_{i},b_{i}|\theta,F_{n})=S_{s}(T_{i},\delta_{i}|\theta,F_{n})+S_{l}(y_{i}|\theta)+S_{r}(b_{i}|\theta),$
(3.15)

where

$S_{l}(y_{i}|\theta)=\frac{\partial}{\partial\theta}E_{b}[\log
p(y_{i}|b_{i};\theta)],$ (3.16)

and

$S_{r}(b_{i}|\theta)=\frac{\partial}{\partial\theta}E_{b}[\log
p(b_{i};\theta)],$ (3.17)

are the score functions for longitudinal and random effect components
respectively. Moreover, $S_{s}(T_{i},\delta_{i}|\theta,F_{n})$ is score
function for the survival component which can be expressed as

$\displaystyle S_{s}(T_{i},\delta_{i}|\theta,F_{n})={}$
$\displaystyle\frac{\partial}{\partial\theta}E_{b}\big{[}\log
p(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F_{n}})\big{]}$ (3.18)
$\displaystyle={}$ $\displaystyle
E_{b}\left[\delta_{i}\left\\{\left(\begin{array}[]{c}m_{i}(T_{i})\\\ \alpha
x_{i}(T_{i})\\\
w_{i}\end{array}\right)-\frac{L_{1,\theta,F_{n}}(T_{i})}{L_{0,\theta,F_{n}}(T_{i})}\right\\}\right]$
$\displaystyle-
E_{b}\left[\int_{0}^{T_{i}}\left\\{\left(\begin{array}[]{c}m_{i}(s)\\\ \alpha
x_{i}(s)\\\
w_{i}\end{array}\right)-\frac{L_{1,\theta,F_{n}}(s)}{L_{0,\theta,F_{n}}(s)}\right\\}\exp\left\\{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(s)\right\\}d\hat{\Lambda}_{\theta,F_{n}}(s)\right].$

where

$\displaystyle L_{0,\theta,F_{n}}(u)$ $\displaystyle=$ $\displaystyle
E_{F_{n}}\left[Y(u)E_{b}\exp\left\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\right\\}\right],$
$\displaystyle L_{1,\theta,F_{n}}(u)$ $\displaystyle=$ $\displaystyle
E_{F_{n}}\left[Y(u)E_{b}\left(\begin{array}[]{c}m(u)\\\ \alpha x(u)\\\
w\end{array}\right)\exp\left\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\right\\}\right].$
(3.22)

Now we will calculate the score function $B(T_{i},\delta_{i}|\theta,F)$, which
is Hadamard differentiable with respect to $F$. For an integrable function $h$
with the same domain as $F$, we can express

$\displaystyle B(T_{i},\delta_{i}|\theta,F)h={}$ $\displaystyle
d_{F}E_{b}\big{[}\log
p(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F})\big{]}h$
$\displaystyle={}$ $\displaystyle
E_{b}\bigg{\\{}\delta_{i}\bigg{[}\frac{E_{h}dN(T_{i})}{E_{F}dN(T_{i})}-\frac{E_{h}Y(T_{i})E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(T_{i})\\}]}{E_{F}Y(T_{i})E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(T_{i})\\}]}\bigg{]}$
$\displaystyle-\exp\\{{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(T_{i})}\\}\int_{0}^{T_{i}}\frac{E_{h}dN(u)}{E_{F}Y(u)E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\\}]}$
$\displaystyle+\exp\\{{\gamma^{\prime}w_{i}+\alpha^{\prime}m_{i}(T_{i})}\\}\int_{0}^{T_{i}}\frac{E_{F}dN(u)E_{h}Y(u)E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\\}]}{\big{(}E_{F}Y(u)E_{b}[\exp\\{\gamma^{\prime}w+\alpha^{\prime}m(u)\\}]\big{)}^{2}}\bigg{\\}}.$

where, $d_{F}\log P(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F})$
represents the Hadamard derivative of $\log
P(T_{i},\delta_{i}|b_{i};\theta,\hat{\Lambda}_{\theta,F})$ with respect to $F$
(hirose2011efficiency).

Theorem 1: At the true value of $(\theta,F)$, we are going to proof the
followings

1\. $\hat{\Lambda}_{\theta_{0},F_{0}}(t)=\Lambda_{0}(t)$, the true cumulative
hazard and

2\. The score function $S(T,\delta,y,b|\theta_{0},F_{0})$ defined in (3.15) is
the efficient score function where we drop the subscript $i$.

Proof: Replace $F_{n}$ by $F_{0}$,we get from (3.10)

$\hat{\Lambda}_{\theta_{0},F_{0}}(t)=\int_{0}^{t}\frac{E[dN(u)]}{E[Y(u)E_{b}\exp\\{\gamma_{0}^{\prime}w+\alpha_{0}^{\prime}m(u)\\}]}.$
(3.23)

where $E$ is the expectation with respect to the true distribution $F_{0}$. At
the true value of the parameters $(\theta,F)$ we can write

$E[dN(u)]=E[Y(u)E_{b}\exp\\{\gamma_{0}^{\prime}w+\alpha_{0}^{\prime}m(u)\\}].$
(3.24)

So from this point of view, we have
$\hat{\Lambda}_{\theta_{0},F_{0}}(t)=\Lambda_{0}(t).$

The score function
$S(T,\delta,y,b|\theta,F)=S_{s}(T,\delta|\theta,F)+S_{l}(y|\theta)+S_{r}(b|\theta)$
in (3.15) has three parts. We know that the longitudinal and random effect
parts are parametric models those do not involve $\Lambda$, so we will work on
the survival part of the score function. So the score function for the
survival part at the true value of the parameters $(\theta,F)$ can be
expressed as

$\displaystyle S^{*}_{\theta}=S_{s}(T,\delta|\theta_{0},F_{0})={}$
$\displaystyle E_{b}\left[\delta\left\\{\left(\begin{array}[]{c}m(T)\\\
\alpha_{0}x(T)\\\
w\end{array}\right)-\frac{L_{1}(T)}{L_{0}(T)}\right\\}\right]$ (3.25)
$\displaystyle-E_{b}\left[\int_{0}^{T}\left\\{\left(\begin{array}[]{c}m(s)\\\
\alpha_{0}x(s)\\\
w\end{array}\right)-\frac{L_{1}(s)}{L_{0}(s)}\right\\}\exp\left\\{\gamma_{0}^{\prime}w+\alpha_{0}^{\prime}m(s)\right\\}d\Lambda_{0}(s)\right].$

where

$\displaystyle L_{0}(u)$ $\displaystyle=$ $\displaystyle
E\left[Y(u)E_{b}\exp\left\\{\gamma_{0}^{\prime}w+\alpha_{0}^{\prime}m(u)\right\\}\right],$
$\displaystyle L_{1}(u)$ $\displaystyle=$ $\displaystyle
E\left[Y(u)E_{b}\left(\begin{array}[]{c}m(u)\\\ \alpha_{0}x(u)\\\
w\end{array}\right)\exp\left\\{\gamma_{0}^{\prime}w+\alpha_{0}^{\prime}m(u)\right\\}\right].$
(3.29)

which is the efficient score function based on profile-likelihood approach.
The calculation of the efficient score function based on the projection theory
is given in Supplementary Materials (equation *).

###  3.3 Asymptotic Normality of the MLE

Assumptions:

To show the asymptotic normality of the MLE and its asymptotic variance, we
have to consider some assumptions. On the set of cdf functions $\digamma$, we
use the sup-norm, i.e., for $F,F_{0}\in\digamma$,

$||F-F_{0}||_{\infty}=\sup_{u}|F(u)-F_{0}(u)|.$

For $\rho>0$, let

$\zeta_{\rho}=\\{F\in\digamma:||F-F_{0}||_{\infty}<\rho\\}.$

The assumptions are given below

A1: $P(T,\delta|b;\theta_{0},F_{0})>\delta>0$ for some positive constant
$\delta>0$ and $S(\tau)=P(T>\tau)=E[Y(\tau)]>0$.

A2: The range of $x$, $z$ and $w$ are bounded and
$\theta=(\alpha,\beta,\gamma)$ is in the compact set $\Theta$ which follows
$||x||\leq M$, $||z||\leq M$, $||w||\leq M$ and $||\beta||\leq M$,
$||\alpha||\leq M$, $||\gamma||\leq M$ for some $0<M<\infty$.

A3: The empirical cdf $F_{n}$ is $\sqrt{n}$ consistent i.e.
$\sqrt{n}|F_{n}-F_{0}|=O_{p}(1)$.

A4: The efficient information matrix
$I_{s}^{*}=E[S^{*}_{\theta}S^{*^{\prime}}_{\theta}]$ is invertible.

Theorem-2: If the assumptions (A1-A4) hold, then

1\. $\hat{\theta}_{n}\stackrel{{\scriptstyle P}}{{\rightarrow}}\theta_{0}$ as
$n\rightarrow\infty$ and 2\.
$\hat{\Lambda}_{\hat{\theta}_{n},F_{n}}-\Lambda_{0}=o_{p}(1)$.

The proof of Theorem-2 is given in Supplementary Materials.

Theorem-3: The score functions $S_{s}(T,\delta|\theta,F)$ and
$B(T,\delta|\theta,F)$ are defined previously. Suppose for $(A1)$-$(A4)$,
$\hat{\theta}_{n}\stackrel{{\scriptstyle P}}{{\rightarrow}}\theta_{0}$ and
$F_{n}\stackrel{{\scriptstyle P}}{{\rightarrow}}F_{0}$ as
$n\rightarrow\infty$, then we have

$E\bigg{[}\sqrt{n}\bigg{\\{}S_{s}(T,\delta|\hat{\theta}_{n},F_{0})-S_{s}(T,\delta|\theta_{0},F_{0})\bigg{\\}}\bigg{]}=-E\bigg{[}S_{s}(T,\delta|\theta_{0},F_{0})S^{\prime}_{s}(T,\delta|\theta_{0},F_{0})\bigg{]}\bigg{\\{}\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\bigg{\\}}+o_{p}(1),$

and

$\displaystyle
E\bigg{[}\sqrt{n}\bigg{\\{}S_{s}(T,\delta|\hat{\theta}_{n},F_{n})-S_{s}(T,\delta|\hat{\theta}_{n},F_{0})\bigg{\\}}\bigg{]}={}$
$\displaystyle-E\bigg{[}S_{s}(T,\delta|\theta_{0},F_{0})B(T,\delta|\theta_{0},F_{0})\bigg{]}\bigg{\\{}\sqrt{n}(F_{n}-F_{0})\bigg{\\}}$
$\displaystyle+o_{p}\big{(}1+\sqrt{n}(\hat{\theta}_{n}-\theta_{0})\big{)}.$

Remark: The results are obtained without assuming the derivative of the score
functions $\frac{\partial}{\partial\theta}S_{s}(T,\delta|\theta,F)$ and
$d_{F}B(T,\delta|\theta,F)$ exist. This result give us asymptotic expansion of
profile likelihood without differentiating the score function that involve
implicit function.

The proof of Theorem-3 is given in Supplementary Materials.

Theorem-4: If the assumptions $\\{A1,A2,A3,A4\\}$ are satisfied, then a
consistent estimator $\hat{\theta}_{n}$ to the estimating equation

$\sum_{i=1}^{n}S_{s}(T_{i},\delta_{i}|\hat{\theta}_{n},F_{n})=0,$

is an asymptotically linear estimator for $\theta_{0}$ (Hirose) with the
efficient influence function
$(I_{s}^{*})^{-1}S_{s}(T,\delta|\theta_{0},F_{0})$, so that

$\sqrt{n}(\hat{\theta}_{n}-\theta_{0})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(I_{s}^{*})^{-1}S_{s}(T_{i},\delta_{i}|\theta_{0},F_{0})+o_{p}(1)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}N\\{0,(I_{s}^{*})^{-1}\\},$

where $N\\{0,(I_{s}^{*})^{-1}\\}$ is a normal distribution with mean zero and
variance $(I_{s}^{*})^{-1}$. So the estimator $\hat{\beta}_{n}$ is efficient.

The proof of Theorem-4 is given in Supplementary Materials.

##  4 Application

##  5 Discussion

We have proposed a profile likelihood estimation approach for the joint model
and calculated the efficient score function (Theorem-1). Moreover, we showed
the consistency (Theorem-2), asymptotic normality and efficiency of the
profile likelihood estimator (Theorem-3 and Theorem-4). By using the profile
likelihood approach, we got the explicit form of the variance estimator for
profile likelihood estimator. For the numerical results, we have used the aids
data from JM R-package and calculated the standard errors of the estimated
parameters from the efficient score function.

## References

  * Bishop (2006) Bishop, C.M. (2006). Pattern recognition and Machine learning, Springer.
  * Dempster, A., Laird, N., & Rubin, D. (1977) Dempster, A., Laird, N., & Rubin, D.  (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 39, 1–38.
  * Hirose (2011) Hirose, Y. (2011). Efficiency of profile likelihood in semi-parametric models, Ann. Inst. Statist. Math. 63 1247–1275.
  * Hirose (2016) Hirose, Y. (2016). On differentiability of implicitly defined function in semi-parametric profile likelihood estimation, Bernoulli 22 589–614.
  * Hsieh, Tseng & Wang (2006) Hsieh, F. & Tseng, Y.K. & Wang, J.L. (2006). Joint modeling of survival and longitudinal data: likelihood approach revisited. Biometrics 62, 1037–1043.
  * Kalbfleisch & Prentice (2002) Kalbfleisch, J. D. & Prentice, R. L. (2002). The Statistical Analysis of Failure Time Data, 2nd edition, John Wiley & Sons, New York.
  * Kolmogorov et al. (1975) Kolmogorov et al. (1975). Introductory real analysis, Dover, New York.
  * McLachlan & Peel (2000) McLachlan, G. & Peel, D. (2000). Finite Mixture Models, Wiley, New York.
  * McLachlan & Krishnan (2008) McLachlan, G. & Krishnan, T. (2008). The EM Algorithm and Extensions Second Edition, Wiley, New York.
  * Murphy and van der Vaart (2000) Murphy, S.A. and van der Vaart, A.W. (2000) On profile likelihood (with discussion). J. Amer. Statist. Assoc. 95 449–485.
  * Preedalikit et al. (2016) Preedalikit et al. (2016). Joint modeling of survival and longitudinal ordered data using a semiparametric approach. Australian & New Zealand Journal of Statistics  58, 153–172.
  * Rizopoulos (2012) Rizopoulos, D. (2012). Joint Models for Longitudinal and Time-to-Event Data With Applications in R, CRC Press.
  * Seber & Lee (2003) Seber, G.A.F. & Lee, A.J. (2003). Linear Regression Analysis, Second Edition. Wiley, New York.
  * Song, Davidian & Tsiatis (2002) Song, X. & Davidian, M. & Tsiatis, A. A. (2002). A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data. Biometrika 58, 742–753.
  * Tsiatis & Davidian (2001) Tsiatis, A. A. & Davidian, M. (2001). A semiparametric estimator for the proportional hazards model with longitudinal covariates measured with error. Biometrika 88, 447–458.
  * Tsiatis & Davidian (2004) Tsiatis, A. A. & Davidian, M. (2004). Joint modeling of longitudinal and time to event data: an overview. Statistica Sinica 14, 809–834.
  * van der Vaart (1998) van der Vaart, A. W. (1998). Asymptotic Statistics, Cambridge university press, Cambridge.
  * van der Vaart & Wellner (1996) van der Vaart, A. W. & Wellner, J.A. (1996). Weak convergence and Empirical processes, Springer, New York.
  * Wulfsohn & Tsiatis (1997) Wulfsohn, M. S. & Tsiatis, A. A. (1997). A joint model for survival and longitudinal data measured with error. Biometrics 53, 330–339.
  * Zeng & Cai (2005) Zeng, D. & Cai, J. (2005). Asymptotic results for maximum likelihood estimators in joint analysis of repeated measurements and survival time. The Annals of Statistics 33, 2132–2163.