# Synthetic Boost: Leveraging Synthetic Data for Enhanced Vision-Language
Segmentation in Echocardiography

Rabin Adhikari Equal Contribution. The order is based on the result of a rock-
paper-scissors round-robin tournament. 0000-0001-5019-2205 Manish
Dhakal††footnotemark:  0000-0002-0101-5592 Safal Thapaliya††footnotemark:
0000-0002-4463-6700 Kanchan Poudel  0009-0007-7249-5494 Prasiddha Bhandari
0009-0004-0836-8253 Bishesh Khanal  0000-0002-2775-4748

###### Abstract

Accurate segmentation is essential for echocardiography-based assessment of
cardiovascular diseases (CVDs). However, the variability among sonographers
and the inherent challenges of ultrasound images hinder precise segmentation.
By leveraging the joint representation of image and text modalities, Vision-
Language Segmentation Models (VLSMs) can incorporate rich contextual
information, potentially aiding in accurate and explainable segmentation.
However, the lack of readily available data in echocardiography hampers the
training of VLSMs. In this study, we explore using synthetic datasets from
Semantic Diffusion Models (SDMs) to enhance VLSMs for echocardiography
segmentation. We evaluate results for two popular VLSMs (CLIPSeg and CRIS)
using seven different kinds of language prompts derived from several
attributes, automatically extracted from echocardiography images, segmentation
masks, and their metadata. Our results show improved metrics and faster
convergence when pretraining VLSMs on SDM-generated synthetic images before
finetuning on real images. The code, configs, and prompts are available at
https://github.com/naamiinepal/synthetic-boost.

###### Keywords:

Vision-Language Models Vision-Language Segmentation Models Echocardiography
Synthetic Data

##  1 Introduction

Echocardiography (heart ultrasound) is an integral diagnostic tool for several
cardiovascular diseases (CVDs). It is widely used because it is cheap,
portable, has no harmful radiation, and has a high temporal resolution (the
ability to see high-definition images in real-time). Accurately estimating
clinically relevant quantitative measures in echocardiography images, such as
cardiac substructure volumes and Ejection Fraction (EF), requires reliable
segmentation algorithms. However, segmenting various parts of the heart is
challenging as the same standard plane image can have diverse appearances
depending on the operator, and the presence of shadows, speckles, strong
attenuation, and low contrast difference among areas of interest in ultrasound
images [1]. Different CNN- and ViT-based [3] U-Net-like models [2, 6, 9, 21]
are the state-of-the-art segmentation models that rely on supervised training
with a relatively large set of annotated echocardiography images. These
segmentation models, however, must be trained on predefined classes that
necessitate retraining or architecture changes (in the final layer) when new
classes are required. It is also challenging to manually intervene in or
inject specific conditioning and make them explicitly benefit from the
spatiotemporal relationships of different foreground structures. Besides, they
lack explainability and are not resilient to distribution shifts.

Recently, Vision-Language Models (VLMs) have been proposed that learn a joint
representation of image and language [4, 8, 10, 13, 19, 22, 29]. VLMs extract
rich supplementary information via image and language prompt pairs,
potentially aiding deep learning models to benefit from the richer
information. VLMs have one encoder each for image and language inputs, and the
encoders are trained together to optimize a joint representation using losses
such as contrastive loss. Vision-Language Segmentation Models (VLSMs) are
adapted from VLMs where a decoder is added and trained on top of pretrained
VLMs to segment the input image while leveraging information provided by
language prompts [16, 20, 26]. However, almost all VLMs are trained using a
large set of natural images, and no VLSMs are trained on an extensive
collection of ultrasound datasets. Although some recent methods show that VLMs
and VLSMs could be finetuned on limited medical data [18], the performance of
these VLSMs is still below the supervised segmentation networks trained and
optimized for specific datasets and foreground masks.

One major challenge to improving VLSMs for ultrasound images is the lack of
large language-image paired datasets. To address the limited data problem,
generative models like GANs [5] and diffusion models [7] could generate images
with a distribution closer to the real-world samples. Stojanovski et al. [23]
trained Semantic Diffusion Models (SDMs) [25] on the CAMUS dataset [12] to
generate synthetic cardiac ultrasound images and showed that the segmentation
model trained exclusively on a generated dataset results in a test dice score
of $89.0\pm 2.5$ in the CAMUS dataset. The use of synthetic images has not
been explored for VLSMs. In this work, we explore whether the synthetic images
from SDMs can improve the performance of VLSMs in echocardiography images.

Our primary contributions are as follows.

  1. 1.

We show that the VLSMs, pretrained on natural images, generalize to the real
dataset (CAMUS) when finetuned on SDM-generated echocardiography images.

  2. 2.

We show that although numerous synthetic samples alone are not as good as a
small number of real annotated data, the model finetuned on synthetic data is
a good starting point for VLSMs to further finetune on real datasets.

##  2 Methodology

###  2.1 Vision-Language Segmentation Models (VLSMs)

Figure 1:  The basic architecture of CRIS and CLIPSeg VLSMs. The key
components in the architecture are a Text Encoder, an Image Encoder, a Vision-
Language Decoder (VLD), and an Aggregator. The images and the corresponding
prompts are passed to the CLIP image and text encoders, respectively. The
Aggregator generates intermediate representations utilizing image-level,
sentence-level, or word-level representations to feed to the VLD. The VLD
outputs a binary segmentation mask for an image-text pair.

CLIP [19] is a widely used VLM that jointly trains an image encoder and a text
encoder to project semantically similar image-text pairs closer together and
semantically disjoint image-text pairs farther apart. As shown in Fig. 1, the
contrastive feature representation obtained from the two encoders of CLIP is
fed to a vision-language decoder which generates a binary segmentation mask.
We investigate CLIPSeg [16] and CRIS [26], two state-of-the-art VLSMs for
natural images, with various combinations of language prompts, real images,
and synthetic images.

We use the publicly accessible CLIPSeg and CRIS weights learned during
pretraining on natural image-text pairings [11, 28]. The two VLSMs are
finetuned on echocardiography datasets, starting from their publicly available
pretrained weights. The two echocardiography datasets are: (i) CAMUS [12], and
(ii) SDM CAMUS [23]. To test if pretraining with synthetic data boosts the
segmentation performance on the real data, natural images pretrained VLSMs are
further trained on extensive synthetic data and then finetuned on a smaller
CAMUS dataset.

###  2.2 Datasets

####  2.2.1 CAMUS

CAMUS [12] is a cardiac segmentation dataset containing 2D apical two-chamber
(2C) and four-chamber (4C) views from $500$ patients at both end-diastole (ED)
and end-systole (ES) cycles. The dataset contains the semantic segmentation of
the left ventricular cavity, the myocardium, and the left atrial cavity. The
original dataset randomly sampled images from $50$ patients as the official
test split, and the remaining $450$ kept in the train split. From those
remaining $450$ patients, like Stojanovski et al. [23], we selected the first
$50$ patients for validation and the remaining $400$ for the training. The
number of train/val/test images is $1,600$/$400$/$200$.

####  2.2.2 Synthetic Echocardiography

We use the synthetic echocardiography images proposed by Stojanovski et al.
[23], generated using SDMs [25]. This model takes perturbed anatomical masks
as conditioning information to denoise the noisy images and generates
echocardiographic images. Our experiments use $9,000$ synthetic images
($8,000$ for training and $1,000$ for validation) provided by the authors, the
same splits they used to train and validate a vanilla U-Net [21] model.

###  2.3 Prompt Engineering

####  2.3.1 Prompts for CAMUS

For our experiments, prompts, images, and masks are needed as triplets but are
unavailable in the CAMUS dataset. Finding the best prompt for the task is
challenging, and creating the prompts manually for each image and mask pair is
tedious and not scalable when the dataset size increases. Also, the choice of
prompts seems to significantly affect the performance of the VLMs in the
medical domain [18].

Table 1: The description of the attribute and its possible values. The prompt number aside shows the prompt in which the attribute is introduced. | Description | Possible Values
---|---|---
P0 | Empty String |
P1 | Target Structure | left ventricular cavity, myocardium, or left atrium cavity
P2 | Apical View | two-chamber view or four-chamber view
P3 | Cardiac Cycle | end of systole or diastole cycle
P4 | Patient’s Sex | male or female
P5 | Patient’s Age | all ages
P6 | Image Quality | good, medium, or poor
P7 | Structure’s Shape | circle, triangle, oval, square, or rectangle

We follow Poudel et al. [17] to generate automatic prompts adapted for the
CAMUS dataset to explore if specific image features could be aligned to
language prompts explaining those features. The foreground cardiac structure’s
size and shape depend on the subjects’ age, sex, and cardiac cycle phase.
Similarly, image quality information may help models adapt accordingly. As
shown in Table 1, various language prompts are designed by including words
corresponding to the target structure name, its shape, the information about
apical views, cardiac cycle phase, the subject’s sex, the subject’s age, and
image quality (labeled by an expert within the CAMUS dataset). There are $7$
attributes generated for the CAMUS dataset and $7$ prompts (P1 - P7) from the
attributes, each added incrementally. P0 is an empty string. The attributes in
P1 - P7 are ordered in descending order of the attribute’s perceived
importance (P1 being the most important).

The sources of the attributes are listed below.

  1. 1.

Image Filename: We parse the images’ filenames and masks to get the anatomical
structure to segment, apical view, and cardiac cycle.

  2. 2.

Image Metadata: We parse the official metadata provided with the images and
masks to get patients’ sex, age, and image quality.

  3. 3.

VQA Model: We use OFA (One For All) VQA [24] to get target structures’ shapes.
The VQA model is presented with the question, What is the shape of the
<structure> in the green box?. Here, the green box is the boundary of the
target structure extracted from its mask.

One example prompt P7 with seven attributes: Left ventricular cavity of oval
shape in two-chamber view in the cardiac ultrasound at the end of the diastole
cycle of a 40-year-old female with poor image quality.

####  2.3.2 Prompts for SDM CAMUS

We did not use the image quality attribute in SDM CAMUS dataset as the
synthetic images’ quality is not annotated. When synthesizing the prompts, we
used the SDM CAMUS dataset’s values derived from the original dataset for all
other attributes: patient id, view information, and cardiac cycle. One example
prompt P6 for the SDM CAMUS dataset: Left ventricular cavity of oval shape in
two-chamber view in the cardiac ultrasound at the end of the diastole cycle of
a 40-year-old female.

##  3 Experimental Settings

Unless specified, the VLSM’s hyperparameters are the same as mentioned in the
original implementation by the respective authors for all experiments. The
models are finetuned and inferred in NVIDIA GeForce RTX 3090, Titan Xp, and
V100 GPUs. We use float-16 mixed-precision training for models with different
batch sizes of $32$ and $128$ for CRIS and CLIPSeg, respectively. The batch
sizes were chosen to utilize the full memory of the GPUs (maximum 24GB); since
CRIS has a greater memory footprint than CLIPSeg, we reduced the former’s
batch size.

We use AdamW [15] optimizer with the weight decay of $10^{-3}$ and an initial
learning rate of $2\times 10^{-3}$ and $2\times 10^{-5}$ for CLIPSeg and CRIS,
respectively. The learning rate is reduced by a factor of $10$ if validation
loss does not decrease for $5$ consecutive
epochs111https://pytorch.org/docs/stable/generated/torch.optim.lr˙scheduler.ReduceLROnPlateau.html.

Three different strategies are employed to train the models: (i) training only
on real data from the CAMUS dataset (real), (ii) training only on synthetic
data generated from the SDM model (synthetic), and (iii) first training the
model on the synthetic data, then finetuning on the real data (synth-PT:real-
FT). CLIPSeg and CRIS resize the input images to $416\times 416$ and
$352\times 352$, respectively. We normalize the resized images with the means
and standard deviations provided by the respective models. No augmentation and
further post-processing are done to assess the models’ raw performance.

We used the weighted sum of soft Dice and Binary Cross Entropy losses with
weights $1$ and $0.2$, respectively. All the dice scores are computed at
$512\times 512$ (nearly the median width of the dataset), resizing the model’s
output when required. For each experiment, the metrics reported are for the
model with the best dice score on the validation set, across the epochs, with
an output threshold of $0.5$ on the predicted binary segmentation map.

To study the ability of the VLSMs to represent the alignment of image-text
pairs, we perform two experiments: (i) freezing the VLM encoders of CRIS and
CLIPSeg, and (ii) unfreezing the VLM encoders during finetuning on all
datasets. The dice score for the unfrozen encoders is shown in the Table 2
whereas that of the frozen ones is demonstrated in Table 3.

##  4 Results

###  4.1 Synthetic data is better than no data

Table 2 shows that while the VLSMs pretrained on natural images perform very
poorly on ultrasound images in zero-shot segmentation, models trained on
synthetic data provide much better results in real ultrasound images.

Table 2: The dice score (mean $\pm$ std) of models trained and validated using
various strategies and evaluated on the CAMUS’s official test split when the
encoders of the VLMs are unfrozen. The zero-shot performance of the models is
extracted from Poudel et al. [17] for comparison.

Strategy | Prompt $\rightarrow$ | P0 | P1 | P2 | P3 | P4 | P5 | P6 | P7
---|---|---|---|---|---|---|---|---|---
Model $\downarrow$
zeroshot | CLIPSeg | $0.00\scriptscriptstyle\pm 0.0$ | $0.00\scriptscriptstyle\pm 0.0$ | $0.21\scriptscriptstyle\pm 1.8$ | $0.16\scriptscriptstyle\pm 1.9$ | $0.19\scriptscriptstyle\pm 2.1$ | $0.51\scriptscriptstyle\pm 3.7$ | $0.46\scriptscriptstyle\pm 3.1$ | $1.81\scriptscriptstyle\pm 6.6$
CRIS | $23.53\scriptscriptstyle\pm 12.0$ | $9.04\scriptscriptstyle\pm 13.9$ | $8.36\scriptscriptstyle\pm 13.2$ | $8.24\scriptscriptstyle\pm 13.2$ | $8.24\scriptscriptstyle\pm 13.2$ | $8.24\scriptscriptstyle\pm 13.2$ | $8.24\scriptscriptstyle\pm 13.2$ | $5.45\scriptscriptstyle\pm 10.4$
synthetic | CLIPSeg | $45.69\scriptscriptstyle\pm 13.2$ | $84.24\scriptscriptstyle\pm 12.0$ | $84.87\scriptscriptstyle\pm 10.9$ | $85.27\scriptscriptstyle\pm 9.7$ | $84.38\scriptscriptstyle\pm 11.0$ | $83.18\scriptscriptstyle\pm 12.8$ | $83.32\scriptscriptstyle\pm 12.5$ | N/A
CRIS | $42.29\scriptscriptstyle\pm 17.6$ | $84.72\scriptscriptstyle\pm 11.9$ | $84.72\scriptscriptstyle\pm 10.5$ | $85.48\scriptscriptstyle\pm 10.2$ | $85.12\scriptscriptstyle\pm 11.2$ | $85.84\scriptscriptstyle\pm 10.0$ | $84.35\scriptscriptstyle\pm 13.3$ | N/A
real | CLIPSeg | $\mathbf{46.52\scriptscriptstyle\pm 13.3}$ | $88.53\scriptscriptstyle\pm 7.2$ | $88.81\scriptscriptstyle\pm 7.2$ | $88.77\scriptscriptstyle\pm 7.2$ | $88.58\scriptscriptstyle\pm 7.7$ | $88.27\scriptscriptstyle\pm 7.4$ | $88.45\scriptscriptstyle\pm 7.5$ | $88.16\scriptscriptstyle\pm 8.0$
CRIS | $46.46\scriptscriptstyle\pm 13.1$ | $91.00\scriptscriptstyle\pm 6.3$ | $91.03\scriptscriptstyle\pm 6.2$ | $89.9\scriptscriptstyle\pm 7.6$ | $90.94\scriptscriptstyle\pm 6.6$ | $90.87\scriptscriptstyle\pm 6.4$ | $90.79\scriptscriptstyle\pm 7.1$ | $90.99\scriptscriptstyle\pm 6.3$
synth-PT: real-FT | CLIPSeg | $46.26\scriptscriptstyle\pm 13.2$ | $88.56\scriptscriptstyle\pm 7.5$ | $89.44\scriptscriptstyle\pm 6.9$ | $89.8\scriptscriptstyle\pm 6.8$ | $88.68\scriptscriptstyle\pm 7.5$ | $88.55\scriptscriptstyle\pm 7.4$ | $89.36\scriptscriptstyle\pm 6.8$ | $89.53\scriptscriptstyle\pm 6.6$
CRIS | $41.09\scriptscriptstyle\pm 18.6$ | $\mathbf{91.26\scriptscriptstyle\pm 6.1}$ | $\mathbf{91.39\scriptscriptstyle\pm 5.9}$ | $\mathbf{91.12\scriptscriptstyle\pm 6.3}$ | $\mathbf{91.04\scriptscriptstyle\pm 7.2}$ | $\mathbf{91.23\scriptscriptstyle\pm 6.4}$ | $\mathbf{91.11\scriptscriptstyle\pm 6.8}$ | $\mathbf{91.08\scriptscriptstyle\pm 6.6}$

###  4.2 Real data is better than synthetic data

Figure 2: Difference in mean dice scores between different training strategies
for CLIPSeg and CRIS for different prompts, relative to real. Pretraining on
synthetic data before finetuning them on real data helps to improve the
performance of VLSMs.

Fig. 2 shows that VLSMs have better dice scores when finetuned in real data
than finetuning only in synthetic data. When comparing the best dice scores
for both strategies, the models trained on the synthetic dataset have a lower
dice score ($-5.19$), which is statistically highly significant by the
Wilcoxon signed-rank test [27] with a p-value of $8.8\times 10^{-73}$, on the
official test split of real images.

###  4.3 Pretraining on synthetic data helps in finetuning on real data

In both CRIS and CLIPSeg, the pretraining on synthetic data and then
finetuning on real data (synth-PT:real-FT strategy) performs better than the
experiments trained with either real or artificial images as illustrated in
Fig. 2. This second stage pretraining strategy has a higher dice score
($+0.34$), which is statistically significant by the Wilcoxon signed-rank test
[27] with a p-value of $8.3\times 10^{-6}$, than the models that haven’t seen
synthetic data.

###  4.4 Unfreezing VLM encoders during finetuning affects models differently

For P0 (empty prompt), the output class is ambiguous for the models. From
Table 3, we can infer that CLIPSeg dealt with this obscurity by predicting a
segmentation map of a union of all the classes, while CRIS chose just noise
(all zeros in the case of the last strategy).

Table 3: The dice score (mean $\pm$ std) on the CAMUS’s official test split
when the encoders of the VLMs are frozen.

Strategy | Prompt $\rightarrow$ | P0 | P1 | P2 | P3 | P4 | P5 | P6 | P7
---|---|---|---|---|---|---|---|---|---
Model $\downarrow$
synthetic | CLIPSeg | $45.71\scriptscriptstyle\pm 13.6$ | $84.08\scriptscriptstyle\pm 11.1$ | $84.01\scriptscriptstyle\pm 10.5$ | $84.48\scriptscriptstyle\pm 11.0$ | $84.02\scriptscriptstyle\pm 11.0$ | $84.47\scriptscriptstyle\pm 10.7$ | $85.47\scriptscriptstyle\pm 9.3$ | N/A
CRIS | $35.13\scriptscriptstyle\pm 19.5$ | $84.19\scriptscriptstyle\pm 13.0$ | $84.02\scriptscriptstyle\pm 12.2$ | $84.62\scriptscriptstyle\pm 11.9$ | $84.94\scriptscriptstyle\pm 11.5$ | $84.23\scriptscriptstyle\pm 12.3$ | $80.70\scriptscriptstyle\pm 17.0$ | N/A
real | CLIPSeg | $\mathbf{46.52\scriptscriptstyle\pm 13.2}$ | $88.81\scriptscriptstyle\pm 7.2$ | $89.04\scriptscriptstyle\pm 7.0$ | $88.65\scriptscriptstyle\pm 7.3$ | $89.05\scriptscriptstyle\pm 7.2$ | $88.54\scriptscriptstyle\pm 7.5$ | $88.61\scriptscriptstyle\pm 7.5$ | $88.54\scriptscriptstyle\pm 7.6$
CRIS | $26.84\scriptscriptstyle\pm 16.2$ | $88.41\scriptscriptstyle\pm 8.7$ | $88.71\scriptscriptstyle\pm 8.6$ | $88.62\scriptscriptstyle\pm 8.8$ | $88.55\scriptscriptstyle\pm 8.7$ | $88.48\scriptscriptstyle\pm 8.6$ | $88.85\scriptscriptstyle\pm 8.4$ | $88.40\scriptscriptstyle\pm 9.8$
synth-PT: real-FT | CLIPSeg | $46.5\scriptscriptstyle\pm 13.3$ | $89.07\scriptscriptstyle\pm 7.1$ | $89.09\scriptscriptstyle\pm 7.1$ | $89.24\scriptscriptstyle\pm 6.7$ | $89.24\scriptscriptstyle\pm 6.9$ | $88.91\scriptscriptstyle\pm 7.2$ | $\mathbf{89.12\scriptscriptstyle\pm 7.0}$ | $89.14\scriptscriptstyle\pm 7.0$
CRIS | $0.04\scriptscriptstyle\pm 0.5$ | $\mathbf{89.21\scriptscriptstyle\pm 7.9}$ | $\mathbf{89.54\scriptscriptstyle\pm 7.4}$ | $\mathbf{89.26\scriptscriptstyle\pm 7.5}$ | $\mathbf{89.41\scriptscriptstyle\pm 7.6}$ | $\mathbf{89.34\scriptscriptstyle\pm 7.8}$ | $89.03\scriptscriptstyle\pm 9.2$ | $\mathbf{89.34\scriptscriptstyle\pm 8.2}$

Fig. 3 shows that CRIS’s performance improves when encoders are not frozen
during finetuning. In contrast, CLIPSeg’s performance degrades when the
encoders are unfrozen for the CAMUS dataset (real one), which seems to have
improved when synthetic data is introduced.

Figure 3: Difference between mean dice scores when the encoders are frozen and
when the encoders are trained for different prompts. CRIS’s model performance
improves when the encoders are trained along with the decoder. In contrast,
CLIPSeg’s performance degrades when encoders are trained.

##  5 Discussion

Although the VLSMs do not improve over the state-of-the-art segmentation
models on the CAMUS dataset (according to the
leaderboard222https://www.creatis.insa-lyon.fr/Challenge/camus/results.html.
Updated January 2023, the maximum mean dice is $94.1$ [14]), it is promising
that they are close. Pretraining with the synthetic samples followed by
finetuning in real samples improves the results compared to finetuning on real
examples without synthetic pretraining. One exciting direction to explore in
the future is to train real and synthetic data together while indicating in
the language prompt whether the sample is real or artificial.

The VLSMs pretrained on natural image-language pairs do not seem to have
captured the language-image relationships common in ultrasound images. Thus,
when finetuning the encoders of VLSMs, the performance improved compared to
freezing the encoders and finetuning only the decoder. CRIS’s performance is
always better when the encoders are finetuned for every strategy, but CLIPSeg
only performs better when the synthetic dataset is introduced. Unfrozen
CLIPSeg performing better when the dataset size is increased may be because,
for CRIS, Wang et al. [26] finetuned the CLIP encoders and the vision-language
decoder for the segmentation task, whereas, Lüddecke et al. [16] froze the
encoders for CLIPSeg. Thus, CLIPSeg’s encoder representation is likely not
well adapted for segmentation as our finetuning of the encoder is limited to
only a few thousand samples.

SDM CAMUS [23] is generated by applying random augmentations to the mask of
the CAMUS dataset. As the dataset was developed by utilizing all the labeled
image-mask pairs in the training set, and the images could not be generated
without the corresponding mask, this questions the “synthetic” portion of the
method (or dataset). This dataset does not solve the medical image
segmentation’s limited paired-data availability problem by generating new
examples. Instead, this is more akin to data augmentation, where the existing
annotated set is augmented with a few thousand transformed pairs by perturbing
existing masks and textures. An important direction in the future would be to
find ways to generate aligned synthetic triplets of language, image, and mask
at scale without annotated image-mask pairs.

##  6 Conclusion

Recent VLSMs trained in large image-language pairs of natural photographs
perform close to the state-of-the-art on the CAMUS echocardiography dataset
when finetuned on the automatically generated prompts. Augmenting training
sets with synthetic images generated from SDM improves VLSMs’ performance.
However, using a relatively large number of synthetic data alone is still
inferior to using a relatively small number of real annotated data. This
suggests that more work is needed in generating better synthetic images whose
distribution is closer to the real data distribution for the echocardiography
dataset. Nevertheless, the synthetic data finetuned model checkpoint seems to
be a good starting point for the segmentation models to finetune on the real
dataset, resulting in improved metrics and faster convergence (on an average,
$4.55$ and $1.71$ times faster for CRIS and CLIPSeg, respectively). While
there is a significant potential for VLSMs for ultrasound image segmentation,
there is a need to develop methods that can generate numerous consistent,
realistic, but synthetic triplets of image, language, and segmentation masks
if one wants to leverage the power of VLSMs.

## References

  * [1] Avola, D., Cinque, L., Fagioli, A., Foresti, G., Mecca, A.: Ultrasound medical imaging techniques: a survey. ACM Computing Surveys (CSUR) 54(3), 1–38 (2021)
  * [2] Deng, K., Meng, Y., Gao, D., Bridge, J., Shen, Y., Lip, G., Zhao, Y., Zheng, Y.: Transbridge: A lightweight transformer for left ventricle segmentation in echocardiography. In: Simplifying Medical Ultrasound: Second International Workshop, ASMUS 2021, Held in Conjunction with MICCAI 2021, Strasbourg, France, September 27, 2021, Proceedings 2. pp. 63–72. Springer (2021)
  * [3] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. In: International Conference on Learning Representations (2020)
  * [4] Fürst, A., Rumetshofer, E., Lehner, J., Tran, V.T., Tang, F., Ramsauer, H., Kreil, D., Kopp, M., Klambauer, G., Bitto, A., et al.: Cloob: Modern hopfield networks with infoloob outperform clip. Advances in neural information processing systems 35, 20450–20468 (2022)
  * [5] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. Advances in neural information processing systems 27, 2672–2680 (2014)
  * [6] Hatamizadeh, A., Tang, Y., Nath, V., Yang, D., Myronenko, A., Landman, B., Roth, H.R., Xu, D.: Unetr: Transformers for 3d medical image segmentation. In: Proceedings of the IEEE/CVF winter conference on applications of computer vision. pp. 574–584 (2022)
  * [7] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. Advances in neural information processing systems 33, 6840–6851 (2020)
  * [8] Huang, Z., Zeng, Z., Liu, B., Fu, D., Fu, J.: Pixel-bert: Aligning image pixels with text by deep multi-modal transformers. arXiv preprint arXiv:2004.00849 (2020)
  * [9] Isensee, F., Jaeger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature methods 18(2), 203–211 (2021)
  * [10] Jia, C., Yang, Y., Xia, Y., Chen, Y.T., Parekh, Z., Pham, H., Le, Q., Sung, Y.H., Li, Z., Duerig, T.: Scaling up visual and vision-language representation learning with noisy text supervision. In: International Conference on Machine Learning. pp. 4904–4916. PMLR (2021)
  * [11] Kazemzadeh, S., Ordonez, V., Matten, M., Berg, T.: Referitgame: Referring to objects in photographs of natural scenes. In: Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP). pp. 787–798 (2014)
  * [12] Leclerc, S., Smistad, E., Pedrosa, J., Østvik, A., Cervenansky, F., Espinosa, F., Espeland, T., Berg, E.A.R., Jodoin, P.M., Grenier, T., et al.: Deep learning for segmentation using an open large-scale dataset in 2d echocardiography. IEEE transactions on medical imaging 38(9), 2198–2210 (2019)
  * [13] Li, Y., Liang, F., Zhao, L., Cui, Y., Ouyang, W., Shao, J., Yu, F., Yan, J.: Supervision exists everywhere: A data efficient contrastive language-image pre-training paradigm. In: International Conference on Learning Representations (2021)
  * [14] Ling, H.J., Garcia, D., Bernard, O.: Reaching intra-observer variability in 2-d echocardiographic image segmentation with a simple u-net architecture. In: IEEE International Ultrasonics Symposium (IUS) (2022)
  * [15] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. In: International Conference on Learning Representations (2018)
  * [16] Lüddecke, T., Ecker, A.: Image segmentation using text and image prompts. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 7086–7096 (2022)
  * [17] Poudel, K., Dhakal, M., Bhandari, P., Adhikari, R., Thapaliya, S., Khanal, B.: Exploring transfer learning in medical image segmentation using vision-language models. arXiv preprint arXiv:2308.07706 (2023)
  * [18] Qin, Z., Yi, H.H., Lao, Q., Li, K.: Medical image understanding with pretrained vision language models: A comprehensive study. In: The Eleventh International Conference on Learning Representations (2022)
  * [19] Radford, A., Kim, J.W., Hallacy, C., Ramesh, A., Goh, G., Agarwal, S., Sastry, G., Askell, A., Mishkin, P., Clark, J., et al.: Learning transferable visual models from natural language supervision. In: International conference on machine learning. pp. 8748–8763. PMLR (2021)
  * [20] Rao, Y., Zhao, W., Chen, G., Tang, Y., Zhu, Z., Huang, G., Zhou, J., Lu, J.: Denseclip: Language-guided dense prediction with context-aware prompting. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 18082–18091 (2022)
  * [21] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18. pp. 234–241. Springer (2015)
  * [22] Singh, A., Hu, R., Goswami, V., Couairon, G., Galuba, W., Rohrbach, M., Kiela, D.: Flava: A foundational language and vision alignment model. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 15638–15650 (2022)
  * [23] Stojanovski, D., Hermida, U., Lamata, P., Beqiri, A., Gomez, A.: Echo from noise: synthetic ultrasound image generation using diffusion models for real image segmentation. arXiv preprint arXiv:2305.05424 (2023)
  * [24] Wang, P., Yang, A., Men, R., Lin, J., Bai, S., Li, Z., Ma, J., Zhou, C., Zhou, J., Yang, H.: Ofa: Unifying architectures, tasks, and modalities through a simple sequence-to-sequence learning framework. In: International Conference on Machine Learning. pp. 23318–23340. PMLR (2022)
  * [25] Wang, W., Bao, J., Zhou, W., Chen, D., Chen, D., Yuan, L., Li, H.: Semantic image synthesis via diffusion models. arXiv preprint arXiv:2207.00050 (2022)
  * [26] Wang, Z., Lu, Y., Li, Q., Tao, X., Guo, Y., Gong, M., Liu, T.: Cris: Clip-driven referring image segmentation. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 11686–11695 (2022)
  * [27] Wilcoxon, F.: Individual comparisons by ranking methods. In: Breakthroughs in Statistics: Methodology and Distribution, pp. 196–202. Springer (1992)
  * [28] Wu, C., Lin, Z., Cohen, S., Bui, T., Maji, S.: Phrasecut: Language-based image segmentation in the wild. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 10216–10225 (2020)
  * [29] Zhai, X., Wang, X., Mustafa, B., Steiner, A., Keysers, D., Kolesnikov, A., Beyer, L.: Lit: Zero-shot transfer with locked-image text tuning. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 18123–18133 (2022)


Martje Rave1, Göran Kauermann1, Martje Rave et al. Department of Statistics,
Faculty of Mathematics, Informatics and Statistics, Ludwig-Maximilians-
Universität München, Germany Martje Rave, Chair of Applied Statistics in
Social Sciences, Economics and Business, Department of Statistics, Faculty of
Mathematics, Informatics and Statistics, Ludwig-Maximilians-Universität
München, Ludwigstr. 33 80539 München Germany<EMAIL_ADDRESS>(+49) 89 2180 2248 (+49) 89 2180-5040 The Skellam Distribution revisited -
Estimating the unobserved incoming and outgoing ICU COVID-19 patients on a
regional level in Germany The Skellam Distribution revisited With the
beginning of the COVID-19 pandemic, we realized the need for comprehensive
data collection and its provision to scientists and experts for proper data
analyses. In Germany, the Robert Koch Institute (RKI) has tried to keep up
with this demand for data on COVID-19, but there were (and still are) relevant
data missing that are needed to understand the whole picture of the pandemic.
In this paper, we take a closer look at the severity of the course of COVID-19
in Germany, for which ideal information would be the number of incoming
patients to ICU units. This information was (and still is) not available.
Instead, the current occupancy of ICU units on the district level was reported
daily. We demonstrate how this information can be used to predict the number
of incoming as well as released COVID-19 patients using a stochastic version
of the Expectation Maximisation algorithm (SEM). This in turn, allows for
estimating the influence of district-specific and age-specific infection rates
as well as further covariates, including spatial effects, on the number of
incoming patients. The paper demonstrates that even if relevant data are not
recorded or provided officially, statistical modelling allows for
reconstructing them. This also includes the quantification of uncertainty
which naturally results from the application of the SEM algorithm. EM; Skellam
distribution, Stochastic EM; Imputation; COVID-19; ICU-patients

##  1 Introduction

Albeit its atrocity, in its aftermath, the COVID-19 pandemic has taught
Germany, among many other countries, the shortcomings of inadequate data
availability in its healthcare system. In fact, in Germany, while intensive
care unit (ICU) occupancy at some point was provided by the DIVI (German
interdisciplinary association of Intensive and emergency medicine), the
important numbers of newly hospitalised patients (incoming) and released
patients (outgoing), either cured or deceased, has (until now) not been
included in the database. While this can be criticised since the most relevant
number, which measures the pressure of the disease - the number of incoming
patients - is not officially released, we show in this paper how to
disentangle incoming and outgoing patients from pure occupancy data using
statistical models. This, in particular, allows us to investigate how
hospitalisations depend on time, age and spatial factors.

We assume that admission to and release of the ICU units follow Poisson
distributions with inhomogeneous intensities. Consequently, the changes in ICU
occupancy result from the difference between incoming and outgoing patients.
This in turn gives the framework of the Skellam distribution, originally
introduced by Skellam (1948). The distribution is easily described as
resulting from the difference of two independent Poisson distributed random
variables. This distributional approach has been used in different settings.
For instance, in sports statistics Karlis and Ntzoufras (2009) apply the
distribution for modelling of the goal difference in football games. In
network analysis, Gan and Kolaczyk (2018) and Schneble and Kauermann (2020)
look at network flows while Koopman et al. (2014) utilise the idea to model
financial trades. Further application areas include image analysis when
comparing intensity differences of pixels, see e.g. Hwang et al. (2007), Hwang
et al. (2011) or Hirakawa et al. (2009). Extensions towards bivariate Skellam
processes are provided e.g. in Genest and Mesfioui (2014), see also Aissaoui
et al. (2017). A general discussion on the Skellam distribution and its
application fields is provided in Tomy and Veena (2022). In this paper, we
provide an application of the Skellam distribution for disentangling incoming
and outgoing patients in ICUs.

The occupancy of ICU units was a central component of the COVID-19 pandemic.
Numerous tools have been developed for forecasting the number of patients who
require ICU admission, see e.g. Grasselli et al. (2020), Goic et al. (2021),
COVID et al. (2020) or Farcomeni et al. (2021) to just mention a few. Our
focus in this paper is not on primarily on prediction but on investigating the
risk of admission and how this depends on the infection rates and further
covariates, including spatial components.

Parameter estimation in the Skellam distribution is cumbersome due to its non-
analytic form of the likelihood function. We refer to Lewis et al. (2016) or
Aissaoui et al. (2017) who pursue moment-based estimation. In this paper, we
show that the stochastic expectation maximisation (SEM) algorithm is a
suitable and numerically stable alternative to available estimation routines.
Originally proposed by Celeux et al. (1996), the stochastic version of the EM
algorithm gained interest in recent years, in particular in mixture models,
see e.g. Noghrehchi et al. (2021). We also refer to Nielsen (2000) for
asymptotic results on the algorithm. The EM algorithm relates the estimation
to a missing data problem, which is easily described. We assume that instead
of the complete data with incoming and outgoing patients, we only observe the
changes in occupancy of ICUs. In other words, the exact number of incoming and
outgoing is missing. Replacing these missing numbers iteratively with
simulated numbers, based on the current estimates of the model, provides the
stochastic version of the ‘E’-step. This, in turn, leads to full data, which
allows for standard maximum likelihood estimation of two Poisson processes -
the M step. The algorithm is easily implemented, and Rubin’s rule, Rubin
(1976), provides inference statements.

A particularly interesting attribute that this approach provides is the
simplification of the initial complexity of the problem. We are able to break
our problem down from a fairly complex distributional assumption, with respect
to deriving an association between the infection rates and the number of
incoming and outgoing patients, to land at essentially two iteratively updated
generalised additive models (GAMs) with simulated responses, each response
simultaneously sampled from a joint distribution, comprised of the product of
two Poisson distributions. This allows us to not only circumvent rather
cumbersome calculations and modifications of the first and second derivative
of the Skellam distribution, as, for example, shown by Schneble and Kauermann
(2022), but also almost effortlessly interpret the association between the
number of incoming and outgoing patients and the infection rate.

The paper is structured as follows. In Section 2, we give a detailed data
description. In Section 3, we elaborate on the model approach to our problem,
while in Section 4, we will provide the results of our model approach. A
simple simulation exercise to validate our findings can be found in Section 5,
and in Section 6, we conclude our paper which also includes a discussion of
the shortcomings of our approach.

##  2 Data Description

The database for our analyses consists of two main components; data on
COVID-19 infections and data on the ICU occupancy of COVID-19 patients. The
infections and the ICU occupancy are collected by the German health care
departments, recorded by the Robert Koch Institute (2021) (RKI), the German
federal government agency and scientific institute responsible for health
reporting and disease control, and published by the RKI and DIVI e.V. (2021),
respectively. We here focus on data during the fourth infection wave in
Germany, i.e. from the $2^{nd}$ of October 2021 until the $17^{th}$ of
November 2021, though the method is readily extendable to other time frames.
We visualise the average infection rates overall districts in Figure 1 (left-
hand side).

The RKI collects and publishes data on infections on a daily basis. Due to
privacy protection, the RKI aggregates the number of COVID-19 patients, ICU
and general admission, over NUTS3 districts, European Commission (2021), and
other demographic groups. These, namely, the age categories; ‘0-4’ year-olds,
‘5-14’ year-olds, ‘15-34’ year-olds, ‘35-59’ year-olds, ‘60-79’ year-olds and
‘80+’ year-olds and the sex; ‘male’, ‘female’ and ‘not disclosed’. The
infections are aggregated over the age groups. The data can be directly
downloaded through the ArcGIS website, Robert Koch Institute (2021). The
infection rates per 100.000 inhabitants are then calculated as a weekly
average for each age group. For each district, the infection rate is averaged
over the seven days immediately preceding the respective observed day change
in ICU occupancy.

The data on ICU occupancy is also collected by the RKI and published by DIVI,
(DIVI e.V., 2021). This data is also on a district level, however, the
occupancy can only be differentiated by the number of beds occupied by
patients infected with COVID-19, by the number of beds occupied by patients
not infected with COVID-19 and the number of empty beds, the sum of which is
the overall ICU capacity in a given district on a given date. We here take the
COVID-19 patients into account and visualise the ICU data for one day in
Figure 1 (right-hand side).

Conveniently, both data sets can also be found in the daily updated GitHub
repository maintained by the RKI, Robert Koch Institute (2023).

Figure 1: A: Summary over all districts of the infection rate per 100.000
inhabitants by age group, ‘45-59’ year-olds, ‘60-79’ year-olds and ‘80+’ year-
olds displayed by date, from the $1^{st}$ of October 2021 until the $18^{th}$
of November 2021, B: The maximum capacity of ICU beds per given district over
the time span from the $1^{st}$ of October 2021 until the $18^{th}$ of
November 2021 by district

##  3 Model

###  3.1 Assumption

Let $Y_{(t,d)}$ be the number of COVID-19 ICU patients in a given district $d$
at day $t$. This is the official number issued by DIVI, described above and
freely accessible from the given sources. We define with $I_{(t,d)}$ the
number of incoming patients in district $d$ at day $t$, which is the number of
newly admitted COVID-19 patients in the ICUs located in district $d$.
Accordingly, we denote with $R_{(t,d)}$ the number of released patients,
meaning that they are discharged, deceased or transferred to a non-ICU unit.
We assume both to come from an inhomogeneous Poisson process such that

$\displaystyle I_{(t,d)}$
$\displaystyle\sim\mbox{Poisson}\left(\lambda^{I}_{(t,d)}\right)$ (1)
$\displaystyle R_{(t,d)}$
$\displaystyle\sim\mbox{Poisson}\left(\lambda^{R}_{(t,d)}\right).$ (2)

The explicit modelling of the intensities $\lambda^{I}_{(t,d)}$ and
$\lambda^{R}_{(t,d)}$ is of primary interest and discussed in depth later in
this section. For now, note that (2) is an approximation, and formally we have
a right censored Poisson distribution with $R_{(t,d)}\leq Y_{(t-1,d)}$ since
no more patients can be released than are currently in the ICU. We can omit
this point, though, since, based on the disease, we know that not all patients
are discharged at a time, so the formal censoring does not play any practical
relevance due to a generally small discharge intensity $\lambda_{(t,d)}$.

With these definitions we can now define the difference $\Delta_{(t,d)}$ in
occupancy of COVID-19 ICU patients per district $d$ and day $t$ to the the
previous day $t-1$ through

$\Delta_{(t,d)}=Y_{(t,d)}-Y_{(t-1,d)}=I_{(t,d)}-R_{(t,d)}.$ (3)

Assuming independence for the number of incoming and outgoing ICU COVID-19
patients together with (1) and (2) leads to a Skellam distribution Skellam
(1948).

$\displaystyle\Delta_{(t,d)}\sim\mbox{Skellam}(\lambda^{I}_{(t,d)},\lambda^{R}_{(t,d)}).$
(4)

Before we derive how to estimate the two intensities in (4) we want to discuss
the suitability of the distributional assumptions. Note that the approach
relies on independence of $I_{(t,d)}$ and $R_{(t,d)}$. This would be violated
if discharges of the ICU in $t$ depend on the number of incoming patients in
$t$. A conceivable scenario where $I_{(t,d)}$ and $R_{(t,d)}$ are dependent
results if the ICUs get to their limit capacity and triage of patients is
inevitable. This situation has not been observed in Germany - over the entire
course of the pandemic - so we can argue that assuming independence between
incoming and outgoing patients is reasonable.

Finally, the intensities $\lambda^{I}_{(t,d)}$ and $\lambda^{R}_{(t,d)}$ are
modelled to depend on a set of covariates denoted by $\mathbf{x}_{(t,d)}$ as
well as previous data. To be specific, we set

$\displaystyle\lambda^{I}_{(t,d)}$
$\displaystyle=\exp\left(\eta_{(t,d)}^{I}+s^{I}(t)+h^{I}(\mbox{longitude}_{d},\mbox{latitude}_{d})\right),$
(5) $\displaystyle\lambda^{R}_{(t,d)}$
$\displaystyle=\exp\Big{(}\eta_{(t,d)}^{R}+s^{R}(t)+h^{R}(\mbox{longitude}_{d},\mbox{latitude}_{d})+\underbrace{\log(\sum^{j=t}_{j=t-56}\omega_{j}\hat{I}_{(j,d)})}_{\mbox{=
offset}}\Big{)},$ (6)

where $\eta_{(t,d)}^{I}$ and $\eta_{(t,d)}^{R}$ are the linear combinations of
the covariates included in the models. Namely, the logged infection rates of
the age groups ‘$35-59$’ year-olds, ‘$60-79$’ year-olds and ‘$80+$’ year-olds,
as well as the weekday, included as a categorical variable, with Friday as its
reference category. $s^{I}(t)$ and $s^{R}(t)$ are smooth functions in time,
and $h^{I}(\mbox{longitude}_{d},\mbox{latitude}_{d})$ and
$h^{R}(\mbox{longitude}_{d},\mbox{latitude}_{d})$ are two-dimensional smooth
functions over the coordinates of the centroids of the respective districts.
Note that $I_{(j,d)}$ is not observed, and we, therefore, replace it with its
simulated value from the ‘E’-step. Moreover, the weights $\omega_{j}$ are
fixed and not estimated but instead obtained from duration time models for
COVID-19 patients in ICU units. We make use of the epidemiological bulletin
published by the RKI in 2020, Tolksdorf et al. (2020a), see Figure 8 in the
appendix. Finally, we impose the standard constraints, i.e. that both
$s^{I}(t)$ and $s^{R}(t)$ as well as the spatial effects
$h^{I}(\mbox{longitude}_{d},\mbox{latitude}_{d})$ and
$h^{R}(\mbox{longitude}_{d},\mbox{latitude}_{d})$ integrate out to one. We
refer to Wood (2017) for more details.

###  3.2 SEM algorithm

Note that estimation of models (5) and (6) would be straightforward if we knew
the data $I_{(t,d)}$ and $R_{(t,d)}$. This suggests making use of an EM
algorithm by treating the estimation as a missing data problem. Note that we
observe $\Delta_{(t,d)}$ from which we can ”calculate” $I_{(t,d)}$ and
$R_{(t,d)}$ . Given the data we have

$\displaystyle I_{(t,d)}=\Delta_{(t,d)}+R_{(t,d)}$ (7)

with the additional constraints that both, $I_{(t,d)}\geq 0$ and
$R_{(t,d)}\geq 0$. Hence, based on the data, we have the joint probability
model for incoming and released ICU patients:

$\displaystyle P(I_{(t,d)}=k,R_{(t,d)}=j|\Delta_{(t,d)}=\delta)$
$\displaystyle\propto\left\\{\begin{tabular}[]{cl}$P(I_{(t,d)}=k)\times
P(R_{(t,d)}=k-\delta)$&\mbox{ for } $j=k-\delta$ \mbox{ and } $k\geq\mbox{
max}(\delta,0)$\\\ $0$&\mbox{ otherwise }\end{tabular}\right.$ (10)

with $P(I_{(t,d)}=k)$ and $P(R_{(t,d)}=k-\delta)$ resulting from the Poisson
model (5) and (6), respectively. While model (10) is a clumsy convolution
model which does not simplify to an analytic form, simulation from the model
is simple by just replacing the infinite pairs $k$ for $I_{(t,d)}$ and
$k-\delta$ for $R_{(t,d)}$ by a set of finite pairs, such that the resulting
cumulative probabilities are approximately equal to one. To be specific, we
have

$\displaystyle P(I_{(t,d)}=k,R_{(t,d)}=j$
$\displaystyle|\Delta_{(t,d)}=\delta,\lambda^{I},\lambda^{R})$
$\displaystyle=\lim_{K\to\infty}\frac{P_{\lambda^{I}}(I_{(t,d)}=k)P_{\lambda^{R}}(R_{(t,d)}=k-\delta)}{\sum_{i=1}^{K}P_{\lambda^{I}}(I_{(t,d)}=i)P_{\lambda^{R}}(R_{(t,d)}=i-\delta)}.$
(11)

We approximate this numerically by assuming that either
$P_{\lambda^{I}}({I_{(}t,d)}=k)$, or $P_{\lambda^{R}}(R_{(t,d)}=k-\delta)$ is
sufficiently close to zero at $k>1000$ making the product of the two
distributions sufficiently close to zero, such that the sum of probabilities
for events $k>1000$ can be ignored. This results in the finite approximation

$\displaystyle P(I_{(t,d)}=k,R_{(t,d)}=j$
$\displaystyle|\Delta_{(t,d)}=\delta,\lambda^{I},\lambda^{R})$
$\displaystyle\approx\frac{P_{\lambda^{I}}(I_{(t,d)}=k)P_{\lambda^{R}}(R_{(t,d)}=k-\delta)}{\sum_{i=1}^{1000}P_{\lambda^{I}}(I_{(t,d)}=i)P_{\lambda^{R}}(R_{(t,d)}=i-\delta)}.$
(12)

Numerically this is easily carried out and allows to simulate data pairs
$(I_{(t,d)}^{*},R_{(t,d)}^{*})$ based on the current estimates of the
intensities using (12) as an approximate version of (10). This provides a
stochastic ‘E’-step and leads to a full data set with (simulated) incoming and
(simulated) released patients for all districts and all time points. With the
resulting (simulated) full data set, we can now directly estimate the
intensities in models (5) and (6), which in turn is conducted in the ‘M’ step.
The ‘M’ step can be carried out by fitting two generalized additive Poisson
models using standard software, see Wood (2017).

Iterating between the two steps gives a stochastic version of the EM
algorithm. Each simulation step provides an estimate, and following the
classical EM algorithm, we can easily see that on average, we increase the
(marginal) likelihood in each step. The outline of which is sketched in
Algorithm 2, see Appendix.

###  3.3 Inference based on SEM

Unlike the EM algorithm, where calculating the variance of the estimates is
not straightforward, and one typically relies on Louis’ formula (Louis, 1982),
the stochastic version allows to take the uncertainty due to the missing data
into account. The derivation shows similarities to Rubin’s formula for
imputation, see Rubin (1976). Let the parameter vector of linear and smooth
functions,
$\hat{\boldsymbol{\beta}}^{(*k)}=(\hat{\boldsymbol{\beta}}^{I(*k)^{T}},\hat{\boldsymbol{\beta}}^{R(*k)^{T}})^{T}$,
be the resulting estimate in the $k^{th}$ step of the SEM algorithm. We assume
$k>k_{0}$, where $k_{0}$ refers to the step when convergence can be
guaranteed. The final estimate results through

$\displaystyle\hat{\boldsymbol{\beta}}=\frac{1}{K-k_{0}}\sum_{k=k_{0}+1}^{K}\hat{\boldsymbol{\beta}}^{(*k)}.$
(13)

The variance is estimated via

$\displaystyle\widehat{\mbox{Var}}(\hat{\boldsymbol{\beta}})=$
$\displaystyle\frac{1}{K-k_{0}}\sum_{k=k_{0}+1}^{K}\widehat{Var}(\hat{\boldsymbol{\beta}}^{(*k)})+\frac{1+(K-k_{0})^{-1}}{(K-k_{0})-1}\sum_{j=k_{0}+1}^{K}(\hat{\boldsymbol{\beta}}^{(*k)}-\hat{\boldsymbol{\beta}})(\hat{\boldsymbol{\beta}}^{(*k)}-\hat{\boldsymbol{\beta}})^{T}$
(14)

where $\widehat{Var}(\hat{\boldsymbol{\beta}}^{(*k)})$ is the variance
estimate in the $k$ iteration step, based on the imputed data set. The latter
directly results through the applied fitting algorithm.

##  4 Results

A great advantage of our approach is that we can directly interpret the
estimated association between the included covariates and the incoming
patients and outgoing patients separately. To do so, we look at covariates
containing information on the infection rates for each of the three age groups
and the weekday effects, which are included as a categorical variable using
reference coding, with Friday as its reference level. The estimated
coefficients and their standard deviation, calculated based on Rubin’s
formula, see Equation 14, are provided in Table 1. We use the last three
hundred runs to determine the coefficient estimates through their median, as
well as their variance through the Equation (14). The estimates over the last
300 runs are shown through line plots in Figures 9 and 10 in the Appendix for
the incoming and outgoing patients, respectively.

First, we look at the association between our covariates and the number of
incoming and outgoing patients, as seen in the middle and right column of the
output table, Table 1. Recall that the weekday effect is included in the model
through a categorical variable, with Friday as its reference category. For the
model estimating the number of incoming patients, keeping respectively all
other variables constant, we can observe that there is an increased number of
incoming patients on other weekdays, compared to Friday, whereas on the
weekend, there is a decreased number of patients, compared to Friday. For
outgoing patients, the behaviour is slightly different. On Monday, Thursday,
Saturday and especially Sunday, fewer patients are released compared to
Friday. Conversely, Tuesday and Wednesday seem slightly increased.

The number of incoming and outgoing patients is positively associated with the
infection rates of all age groups. Notably, the strongest effect exists for
the infection rate of 35 to 59-year-olds. This is interesting, bearing in mind
that 60 to 69-year-olds are the predominant age group in intensive care. We
should, however, not omit that there is strong collinearity between the
infection rates themselves, as seen in Figure 6 in the Appendix.

| Incoming | Outgoing
---|---|---
| Estimates | Std. Dev. | Estimates | Std. Dev.
Intercept | -2.28 | 0.10 | -6.41 | 0.12
Monday Effect | 0.12 | 0.05 | -0.21 | 0.06
Tuesday Effect | 0.14 | 0.05 | 0.03 | 0.06
Wednesday Effect | 0.13 | 0.05 | 0.02 | 0.06
Thursday Effect | 0.14 | 0.05 | -0.10 | 0.06
Saturday Effect | -0.02 | 0.05 | -0.09 | 0.06
Sunday Effect | -0.14 | 0.05 | -0.39 | 0.06
Infection 35-59 yo | 0.24 | 0.05 | 0.28 | 0.06
Infection 60-79 yo | 0.07 | 0.05 | 0.07 | 0.05
Infection 80+ yo | 0.11 | 0.02 | 0.10 | 0.02
Table 1: Estimated coefficients and standard deviations presented on the level
of incoming and outgoing patients. The estimates are the exponential of the
median of the coefficient estimates from the 200th run to the 500th run of the
EM algorithm.

We also included smooth functions to estimate both the spatial-, and the
temporal effects. They are included to pick up on additional spatial and
temporal structural dependencies. Let us first look at the smooth effects in
time, as seen in Figure 2. The averaged smooth function over time for incoming
patients (left-hand side) is generally increasing. Evidently, we can see some
fluctuation and there seems to be a fortnightly rhythm within the overall
trend. Here we observe an increase in the number of incoming patients for the
first seven days, then a decrease in the following seven days, followed by a
subsequent increase, and so forth. In contrast, as shown on the right-hand
side of Figure 2, we see a general decrease in the number of outgoing patients
without a biweekly rhythm.

Figure 2: Estimated smooth functions of all runs, over time, rendered by the
GAMs estimating the number of incoming patients (left hand side) and outgoing
patients (right hand side).  Figure 3: Estimated median smooth functions of
the $200^{th}$ until the $500^{th}$ run over space rendered by the generalized
additive models, estimating the number of incoming patients (left-hand side)
and outgoing patients (right-hand side).

Finally, we look at the spatial effects for the incoming patients, see the
left-hand side of Figure 3, and for the outgoing patients, shown with the
right-hand side of Figure 3. There seems to be an increased level of incoming
patients in Saxony (east Germany) and North Rhine-Westphalia (west Germany)
and a slight increase around the larger cities of Germany (Frankfurt,
Stuttgart and Munich, south and southwest of Germany). We observe a similar
structure in the spatial smooth function in the model estimating the outgoing
patients, except for the the strong increase around Saxony. Overall, we see
clear spatial heterogeneity.

Figure 4: Estimated number of incoming and outgoing patients by date from the
$1^{st}$ of October 2021 until the $18^{th}$ of November 2021, as well as the
total number of COVID-19 patients in the ICUs of Germany.

At last, we visualise in Figure 4 the estimated number of incoming and
outgoing patients, summed up over the entirety of Germany, for the observed
time frame. The left-hand axis scales the number of incoming and outgoing
patients, whereas the right-hand axis scales the number of overall ICU
patients with COVID-19. We see that the model picks up the somewhat constant
occupancy, from the $1^{st}$ of October 2021, until the $17^{th}$ of October
2021, in Germany’s ICUs rather well, where the number of incoming and outgoing
patients are estimated to be similar, if not equal. Thereafter, the number of
ICU patients in the ICU increases, around this time, we also observe a higher
estimated number of incoming patients than outgoing patients. It is not
unusual for patients, especially the critically ill, to stay in the ICU for
more than four weeks, making the divergence in estimation for the number of
incoming patients and outgoing patients entirely plausible.

##  5 Simulation

To investigate the modelling approach we chose a simple version to emulate the
data used above. We use one covariate, randomly drawn from a normal
distribution, whose mean and variance are taken from the observed mean and
variance of the logged infection rates of ‘60-79’ year-olds. We choose this
age group as ‘60-69’ year-olds are the predominant group in the German ICUs
during the fourth wave, as described in Section 2. The coefficients for the
simulation are chosen in a way such that the difference in the simulated
incoming and outgoing patients is somewhat similar to the range of the
difference in the observed incoming and outgoing patients, namely $(-24,20)$
in the observed data.

$\displaystyle X\sim N(1.978,1.397)$ (15)

The incoming and outgoing number of patients are then simulated as shown in
equation 5. The estimated coefficients in twenty runs are shown in Figure 5.

$\displaystyle I$ $\displaystyle\sim
Poi(exp(\beta^{in}_{0}+\beta^{out}_{1}X))$ $\displaystyle R$
$\displaystyle\sim Poi(exp(\beta^{out}_{0}+\beta^{out}_{1}X+log(I_{lag})))$
(16)

We see for the confidence intervals of each of the coefficient estimates of
each of the twenty runs include the real coefficient, except for the “Incoming
Intercept” coefficient in the $12^{th}$ simulated data set. Overall, the
simulation confirms that we are able to uncover incoming and outgoing patients
from pure hospitalisations.

Figure 5: The estimated coefficients for twenty simulated data sets. Algorithm
1 Simulation algorithm

$n=1000+1$ $\triangleright$ Simulate a data set with 1000 entries, we need to
start with $n=1000+1$. A seed is also set.

$X$ $\triangleright$ Draw $X_{c}$ as a vector of $n$ entries from a Normal
distribution with mean $1.978$ and variance $1.397$. This will be the
covariate used to simulate the number of incoming and outgoing patients,
$I^{sim}$ and $R^{sim}$. Let the design matrix, $X$, be a matrix of a $1\times
n$ dimensional one vector column bound to $X_{c}$.

$\vec{\beta}_{In}$ and $\vec{\beta}^{Out}$ $\triangleright$ Set parameter
vector association between responses and covariate. Here,
$\vec{\beta}^{In}=(\beta^{In}_{0},\beta^{In}_{1})=(-2.340,0.8)$ and
$\vec{\beta}^{Out}=(\beta^{Out}_{0},\beta^{Out}_{1})=(0.001,0.1)$

$I_{sim}\sim Poi(exp(\vec{\beta}^{In}X^{T}))$ $\triangleright$ For each entry
in $X$ respectively draw a number of incoming patients from a Poisson
distribution, $Poi(exp(\vec{\beta}^{In}X^{T}))$.

$R_{sim}\sim Poi(exp(\vec{\beta}^{Out}X_{Out}^{T}+log(I^{sim}_{Lag}+0.1)))$
$\triangleright$ Let $X_{Out}$ be the last $n-1$ rows of $X$. Let
$I^{sim}_{Lag}$ be the first $n-1$ entries of $I^{sim}$. Draw $n-1$ entries
from the Poisson distribution,
$Poi(exp(\vec{\beta}^{Out}X_{Out}^{T}+log(I^{sim}_{Lag}+0.1)))$.

Return Simulated data $\triangleright$ Join the last $n-1$ entries of
$I^{sim}$ with $R^{sim}$, calculate and join their difference
$I^{sim}-R^{sim}=\delta^{sim}$. Further, join with the last $n-1$ entries of
$X_{c}$.

##  6 Concluding remarks

Overall, in this application of the SEM, we are not only able to simulate
unobserved data but also estimate the association between the weekday effect
and the infection rates and the number of incoming and outgoing patients in a
simple and intuitive manner. We achieve some insight into the estimated
association between the infection rates and the number of incoming and
outgoing patients. Namely, the driving force of the estimated number of
incoming and outgoing patients seems to be the infection rates of ‘35-59’
year-olds. Although we are not able to validate the predictions against the
actual number of incoming and outgoing ICU patients, our findings seem to be
mostly reasonable. Additionally, the SEM estimate the association of the
simulated number of incoming and outgoing ICU patients and the simulated
covariate, which are obviously observed measures in all but one simulation. In
this situation, the SEM seems to be an appropriate application and allows us
to gain a more complete picture of the COVID-19 pandemic, even when dealing
with incomplete information.

## Acknowledgements

We acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG,
Project KA 1188/13-1).

## Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the
research, authorship and/or publication of this article.

## References

  * Aissaoui et al. (2017) Aissaoui, S. A., Genest, C., and Mesfioui, M. (2017).  A second look at inference for bivariate skellam distributions.  _Stat_ , 6(1), 79–87.
  * Celeux et al. (1996) Celeux, G., Chauveau, D., and Diebolt, J. (1996).  Stochastic versions of the em algorithm: an experimental study in the mixture case.  _Journal of statistical computation and simulation_ , 55(4), 287–314.
  * COVID et al. (2020) COVID, I., Murray, C. J., et al. (2020).  Forecasting covid-19 impact on hospital bed-days, icu-days, ventilator-days and deaths by us state in the next 4 months.  _MedRxiv_.
  * DIVI e.V. (2021) DIVI e.V. (2021).  Daily ICU occupancy data for COVID-19 and non-COVID-19 patients.  Available at https://www.divi.de/register/tagesreport.  Accessed: .
  * European Commission (2021) European Commission (2021).  Eurostat europe nuts maps.  Available at https://ec.europa.eu/eurostat/web/nuts/nuts-maps.  Accessed: .
  * Farcomeni et al. (2021) Farcomeni, A., Maruotti, A., Divino, F., Jona-Lasinio, G., and Lovison, G. (2021).  An ensemble approach to short-term forecast of covid-19 intensive care occupancy in italian regions.  _Biometrical journal_ , 63(3), 503–513.
  * Gan and Kolaczyk (2018) Gan, H. L. and Kolaczyk, E. D. (2018).  Approximation of the difference of two poisson-like counts by skellam.  _Journal of Applied Probability_ , 55(2), 416–430.
  * Genest and Mesfioui (2014) Genest, C. and Mesfioui, M. (2014).  Bivariate extensions of skellam’s distribution.  _Probability in the Engineering and Informational Sciences_ , 28(3), 401–417.
  * Goic et al. (2021) Goic, M., Bozanic-Leal, M. S., Badal, M., and Basso, L. J. (2021).  Covid-19: Short-term forecast of icu beds in times of crisis.  _Plos one_ , 16(1), e0245272.
  * Grasselli et al. (2020) Grasselli, G., Pesenti, A., and Cecconi, M. (2020).  Critical care utilization for the covid-19 outbreak in lombardy, italy: early experience and forecast during an emergency response.  _Jama_ , 323(16), 1545–1546.
  * Hirakawa et al. (2009) Hirakawa, K., Baqai, F., and Wolfe, P. J. (2009).  Wavelet-based poisson rate estimation using the skellam distribution.  In _Computational Imaging VII_ , volume 7246, pages 215–226. SPIE.
  * Hwang et al. (2007) Hwang, Y., Kim, J.-S., and Kweon, I.-S. (2007).  Sensor noise modeling using the skellam distribution: Application to the color edge detection.  In _2007 IEEE conference on computer vision and pattern recognition_ , pages 1–8. IEEE.
  * Hwang et al. (2011) Hwang, Y., Kim, J.-S., and Kweon, I. S. (2011).  Difference-based image noise modeling using skellam distribution.  _IEEE transactions on pattern analysis and machine intelligence_ , 34(7), 1329–1341.
  * Karlis and Ntzoufras (2009) Karlis, D. and Ntzoufras, I. (2009).  Bayesian modelling of football outcomes: using the skellam’s distribution for the goal difference.  _IMA Journal of Management Mathematics_ , 20(2), 133–145.
  * Koopman et al. (2014) Koopman, S. J., Lit, R., and Lucas, A. (2014).  The dynamic skellam model with applications.
  * Lewis et al. (2016) Lewis, J. W., Brown, P. E., Tsagris, M., and Brown, M. P. E. (2016).  Package ‘skellam’.
  * Louis (1982) Louis, T. A. (1982).  Finding the observed information matrix when using the em algorithm.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 44(2), 226–233.  ISSN 00359246.  URL http://www.jstor.org/stable/2345828.
  * Nielsen (2000) Nielsen, S. F. (2000).  The stochastic em algorithm: estimation and asymptotic results.  _Bernoulli_ , pages 457–489.
  * Noghrehchi et al. (2021) Noghrehchi, F., Stoklosa, J., Penev, S., and Warton, D. I. (2021).  Selecting the model for multiple imputation of missing data: Just use an ic!  _Statistics in Medicine_ , 40(10), 2467–2497.
  * Robert Koch Institute (2021) Robert Koch Institute (2021).  Daily COVID-19 cases data.  Available at https://www.arcgis.com/home/item.html?id=f10774f1c63e40168479a1feb6c7ca74.  Accessed: .
  * Robert Koch Institute (2023) Robert Koch Institute (2023).  Daily COVID-19 cases data.  Available at https://github.com/robert-koch-institut.  Accessed: .
  * Rubin (1976) Rubin, D. B. (1976).  Inference and missing data.  _Biometrika_ , 63(3), 581–592.
  * Schneble and Kauermann (2020) Schneble, M. and Kauermann, G. (2020).  Estimation of latent network flows in bike-sharing systems.  _Statistical Modelling_ , page 1471082X20971911.
  * Schneble and Kauermann (2022) Schneble, M. and Kauermann, G. (2022).  Estimation of latent network flows in bike-sharing systems.  _Statistical Modelling_ , 22(4), 349–378.  10.1177/1471082X20971911.  URL https://doi.org/10.1177/1471082X20971911.
  * Skellam (1948) Skellam, J. G. (1948).  A probability distribution derived from the binomial distribution by regarding the probability of success as variable between the sets of trials.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 10(2), 257–261.
  * Tolksdorf et al. (2020a) Tolksdorf, K., Buda, S., Schuler, E., Wieler, L. H., and Haas, W. (2020a).  Eine hoehere letalitaet und lange beatmungsdauer unterscheiden covid-19 von schwer verlaufenden atemwegsinfektionen in grippewellen.  (41), 3–10.  http://dx.doi.org/10.25646/7111.
  * Tolksdorf et al. (2020b) Tolksdorf, K., Buda, S., Schuler, E., Wieler, L. H., and Haas, W. (2020b).  Eine höhere letalität und lange beatmungsdauer unterscheiden covid-19 von schwer verlaufenden atemwegsinfektionen in grippewellen.  (41), 3–10.  http://dx.doi.org/10.25646/7111.
  * Tomy and Veena (2022) Tomy, L. and Veena, G. (2022).  A retrospective study on skellam and related distributions.  _Austrian Journal of Statistics_ , 51(1), 102–111.
  * Wood (2017) Wood, S. N. (2017).  _Generalized additive models: An introduction with R_.  Boca Raton: CRC press.

##  Appendix A Sketch of Algorithm

Algorithm 2 EM algorithm

$\delta$ $\triangleright$ $\delta$ is the observed difference.

$\lambda^{0}_{In}$, $\lambda^{0}_{Out}$ $\triangleright$ Choose starting
values for the Poisson parameters of the number of incoming and outgoing
patients, respectively.

$P(I^{0}_{(t,d)}=k,R^{0}_{(t,d)}=k-\delta|\lambda^{0}_{In},\lambda^{0}_{Out})\sim
Poi_{\lambda^{0}_{In}}(k)Poi_{\lambda^{0}_{Out}}(k-\delta)$ $\triangleright$
Simulate $I^{0}_{(t,d)}$ and $R^{0}_{(t,d)}$ from a joint probability
distribution of two Poisson distributions, given the Poisson parameters
previously chosen.

for 500 iterations- $i\in\left\\{1,\dots,500\right\\}$ do

M Step $\triangleright$ Estimate $\hat{\lambda}^{i}_{In}$ and
$\hat{\lambda}^{i}_{Out}$ by using generalised additive models to estimate the
association between the used covariates and the previously simulated number of
incoming and outgoing patients, $I^{i-1}_{(t,d)}$ and $R^{i-1}_{(t,d)}$.

‘E’-step $\triangleright$ Simulate the number of incoming and outgoing
patients through the joint probability distribution.

$P(I^{i}_{(t,d)}=k,R^{i}_{(t,d)}=k-\delta|\lambda^{i}_{in},\lambda^{i}_{out})\sim
Poi_{\hat{\lambda}^{i}_{In}}(k)Poi_{\hat{\lambda}^{i}_{Out}}(k-\delta)$

end for

Return $\triangleright$ A list of estimated parameters (M-Step) and simulated
number of incoming and outgoing patients (‘E’-step) for each iteration.

##  Appendix B Additional Plots

Figure 6: A plot outlining the linear correlation between the infection rates
of ‘35-59’ year-olds, ‘60-79’ year-olds and ‘80+’ year-olds. On the lower
triangle, we see their respective scatter plots, along the diagonal, their
distribution and on the upper triangle, their pairwise Pearson’s correlation
coefficients. Figure 7: Estimated maximum number of incoming A) and outgoing
patients B) per district from the $2^{nd}$ of October 2021 until the $17^{th}$
of November 2021. Figure 8: Percentage of outgoing ICU patients after the day
of admission, as presented by Tolksdorf et al. (2020b). Figure 9: Coefficients
estimated by the generalised additive models of the last three hundred runs of
the EM algorithm of the incoming patients. Figure 10: Coefficients estimated
by the generalised additive models of the last three hundred runs of the EM
algorithm of the outgoing patients.