2. De tijdreeksmethode

2.1 Het model

Voor de berekening van de kwartaalcijfers gaan we uit van het bivariate structurele tijdreeksmodel gegeven door 
\begin{eqnarray*}
\begin{pmatrix} \hat{y}_t^{\rm R} \\ \hat{y}_t^{\rm I} \end{pmatrix} =
\begin{pmatrix} \theta_t \\ \theta_t \end{pmatrix} +
\begin{pmatrix} 0 \\ \lambda_t \end{pmatrix} +
\begin{pmatrix} e_t^{\rm R}  \\ e_t^{\rm I} \end{pmatrix} 
\mbox{met $\theta_t = L_t + S_t$ en kwartaal $t$.}
\end{eqnarray*}

Hierin is 

• \(\hat{y}_t^{\rm R}\): de directe schatting in kwartaal \(t\) gebaseerd op de volledige respons (reguliere reeks).
• \(\hat{y}_t^{\rm I}\): de directe schatting in kwartaal \(t\) gebaseerd op de internetrespons (internetreeks).
• \(\theta_t\): de onbekende populatieparameter in kwartaal \(t\).
• \(L_t\): de trend, gemodelleerd als smooth-trendmodel.
• \(S_t\): de seizoenscomponent met kwartaal als periode, gemodelleerd met een trigonometrisch seizoen model. 
• \(\lambda_t\): het systematisch verschil tussen de internetreeks en de reguliere reeks, gemodelleerd als random walk.
• \(e_t^j\): de meetfout van \(\hat{y}_t^j\) voor \(j \in \{\rm R, \rm I\} \), gemodelleerd als \(e_t^j = \sqrt{\widehat{\rm V}(\hat{y}_t^j) } \tilde{e}_t^j \) waarbij \(\widehat{\rm V}(\hat{y}_t^j)\) de geschatte variantie van \(\hat{y}_t^j\) is en \(\tilde{e}_t^j\) witte ruis.
De volledige uitwerking van dit model staat in Bijlage I.

2.2 De directe schattingen

De directe schattingen \(\hat{y}_t^{\rm R}\) en \(\hat{y}_t^{\rm I}\) kunnen berekend worden vanaf het eerste kwartaal van 2014. Voor de jaren 2014 t/m 2019 zijn deze schattingen gebaseerd op de gewogen jaarbestanden van de gezondheidsenquête. De weging zoals die vanaf 2014 plaatsvindt, inclusief een update met de veranderingen per 2018 is beschreven in Boonstra (2019). 

In 2020 is er voor de berekening van de kwartaalcijfers steeds een gewogen bestand samengesteld door de tot dan toe beschikbare respons via de gegeneraliseerde regressieschatter (GREG, Särndal e.a. 1992) op te hogen naar de populatie van 2020. Hierbij is gebruik gemaakt van het weegmodel uit eerdere jaren, waarbij geen rekening gehouden is met de in 2018 ingevoerde doelgroepenbenadering. De eerste twee kwartalen over 2020 zijn in het derde kwartaal van 2020 berekend en gepubliceerd in augustus. De directe schattingen voor deze twee kwartalen zijn daarom gebaseerd zijn op een gewogen halfjaarsbestand. De directe schattingen voor het derde kwartaal zijn gebaseerd op een gewogen negenmaandsbestand en de directe schattingen voor het vierde kwartaal op het gewogen jaarbestand van 2020. 

De directe schatting \(\hat{y}_t^{\rm R}\) is berekend als gewogen gemiddelde van de volledige respons (cawi plus capi) in kwartaal \(t\) van het beschikbare gewogen responsbestand en de directe schatting \(\hat{y}_t^{\rm I}\) als gewogen gemiddelde van alleen de internetrespons in kwartaal \(t\) van hetzelfde responsbestand. De standaardfouten van de directe schattingen, d.w.z. \( \sqrt{\widehat{\rm V}(\hat{y}_t^{\rm R})} \) en \(\sqrt{\widehat{\rm V}(\hat{y}_t^{\rm I})}\), zijn berekend in R (R Core Team, 2015) met het package ‘Survey’ (Lumley, 2014), waarbij rekening gehouden is met het steekproefontwerp van de gezondheidsenquête. Er wordt maandelijks een zelfwegende gestratificeerde tweetrapssteekproef getrokken met gemeente als cluster en corop als stratum in de eerste trap met een minimale clusteromvang gelijk aan één. In de tweede trap worden personen per gemeente getrokken. In eerdere jaren werd deelgemeente als cluster gebruikt. Bij het schatten van de standaardfouten gaan we uit van een gestratificeerd steekproefontwerp met de combinatie maand × provincie als stratum, aangezien de combinatie maand × corop tot strata met te weinig respons leidt. 

De directe schattingen van de internetreeks over de kwartalen 1 en 2 en de directe schattingen voor de reguliere reeks zijn dus berekend op basis van het gewogen halfjaarsbestand. De directe schattingen van beide reeksen over het vierde kwartaal zijn berekend op basis van het gewogen jaarbestand. 

Voor het derde kwartaal is het gewogen negenmaandsbestand gebruikt gebaseerd op de beschikbare respons van januari t/m september. In het derde kwartaal is de capi-respons maar voor twee derde beschikbaar, omdat de capi-waarneming pas in augustus is opgestart. Er is daarom een correctie toegepast bij het berekenen van de directe schatting van de reguliere reeks over kwartaal 3. De directe schatting over kwartaal 3 is berekend als het gewogen gemiddelde van de volledige respons (cawi plus capi) in augustus en september minus 1/3 maal de breuk \(\hat{\lambda}_t\) zoals die door het tijdreeksmodel geschat is voor kwartaal 2. Op de bijbehorende standaardfouten is geen correctie toegepast. 

Op deze manier zijn de twee inputreeksen voor het tijdreeksmodel samengesteld. De reguliere reeks, bestaande uit de schattingen \(\hat{y}_t^{\rm R}\), loopt van het eerste kwartaal 2014 t/m het vierde kwartaal van 2020, waarbij het tweede kwartaal van 2020 ontbreekt. De internetreeks bestaat uit de schattingen \(\hat{y}_t^{\rm I}\) en loopt van het eerste kwartaal 2014 t/m het vierde kwartaal van 2020.

2.3 De modelgebaseerde schattingen

Aan de hand van het tijdreeksmodel uit paragraaf 2.1 kunnen schattingen gemaakt worden voor de populatieparameter \(\theta_t\) vanaf het eerste kwartaal van 2014 t/m het laatste kwartaal van 2020. Hiertoe berekenen we via het Kalman filter (Durbin en Koopman, 2012) de gefilterde schattingen van \(\hat{L}_t + \hat{S}_t\). Gefilterd wil zeggen dat de schatting \(\hat{L}_t + \hat{S}_t\) in kwartaal \(t\) gebaseerd is op alle beschikbare respons vanaf het eerste kwartaal van 2014 t/m kwartaal \(t\). De nowcasts voor de reguliere reeks in het tweede kwartaal, waar capi ontbreekt, volgen hier dan automatisch uit. Hierbij worden twee sterke aannames gemaakt. Ten eerste wordt verondersteld dat de samenstelling van de internetrespons niet verandert tijdens de coronacrisis. Deze aanname is geëvalueerd via een responsanalyse en dit lijkt inderdaad het geval te zijn. De tweede veronderstelling is dat het verschil tussen de capi- en cawi-respons niet verandert door de coronacrisis. Deze modelveronderstelling kan niet worden geëvalueerd.

Om het Kalman filter toe te kunnen passen wordt het model uit paragraaf 2.1 eerst als toestandsruimtemodel geschreven:
\begin{align*}
& \textbf{y}_t = \textbf{Z}_t \alpha_t \\
&  \alpha_t  =  \textbf{T} \alpha_{t-1} + \eta_t \mbox{ voor kwartaal $t$, waarbij} \\
& \eta_t  \sim  {\cal N}(0,\textbf{H}). 
\end{align*}

In het toestandsruimtemodel wordt de vector \(\textbf{y}_t = (\hat{y}_t^{\rm R},\hat{y}_t^{\rm I})’ \) van de directe schattingen uitgedrukt in een vector van de niet-waargenomen toestandsvariabelen \(\alpha_t\), zoals de trend en de seizoencomponenten. Verder wordt in het toestandsruimtemodel het verloop van deze toestandsvariabelen over de tijd beschreven. De covariantiematrix \(\textbf{H}\) van de vector \(\eta_t\) met storingstermen bevat de hyperparameters van het model, die geschat worden met maximum likelihood (ML). De berekeningen zijn uitgevoerd in Ssfpack 3.0 (Koopman et al., 2008) in combinatie met Ox (Doornik, 2009).

Bij het berekenen van de cijfers over de eerste twee kwartalen en bij het berekenen van het derde en vierde kwartaal zijn voor alle doelvariabelen de modelaannames geëvalueerd door voor zowel de reguliere reeks als de internetreeks de gestandaardiseerde innovaties te analyseren. Innovaties zijn de fouten van de voorspellingen van de waarnemingen \(\textbf{y}_t = (\hat{y}_t^{\rm R},\hat{y}_t^{\rm I})’ \) gegeven door het Kalman filter gebaseerd op gegevens die tot en met een kwartaal eerder beschikbaar zijn. De door het Kalman filter voorspelde waarden noteren we met \( \textbf{y}_{t|t-1} = (\hat{y}_{t|t-1}^{\rm R},\hat{y}_{t|t-1}^{\rm I})’ \). De innovaties voor de reguliere reeks zijn dan gegeven door \(\hat{y}_{t}^{\rm R} - \hat{y}_{t|t-1}^{\rm R}\) en voor de internetreeks door \(\hat{y}_{t}^{\rm I} - \hat{y}_{t|t-1}^{\rm I}\). De gestandaardiseerde innovaties zijn gelijk aan de innovaties gedeeld door de bijbehorende standaardfouten. Onder het model zijn de gestandaardiseerde innovaties standaard normaal verdeeld. Via deze gestandaardiseerde innovaties is het model onderzocht op normaliteit, heteroscedasticiteit en autocorrelatie. Hiervoor zijn een aantal toetsen uitgevoerd, zie Durbin en Koopman, (2012) Hoofdstuk 2 voor details. 

Bij het toepassen van het tijdreeksmodel blijken bij enkele variabelen in eerdere jaren uitbijters voor te komen. Het gaat om huisartscontact in kwartaal 3 van 2016, dagelijks roken in kwartaal 2 van 2019 en overgewicht in kwartaal 4 van 2018. De gestandaardiseerde innovaties vallen in deze kwartalen buiten het interval (-2,2). Bij het berekenen van de eerste twee kwartalen worden bij ervaren gezondheid en de drie variabelen over het zorggebruik in de eerste kwartalen van 2020 gestandaardiseerde innovaties gemeten die in absolute waarde variëren van 3 tot 6. Dat is een teken dat voor deze variabelen het tijdreeksmodel de effecten van corona in 2020 in de eerste twee kwartalen niet goed beschrijft en aangepast moet worden. 

2.4 Aanpassing model vanwege corona-effecten

Schattingen gebaseerd op het tijdreeksmodel lenen informatie uit het verleden om zo de nauwkeurigheid van de schattingen te verbeteren. Er wordt daarbij aangenomen dat de cijfers uit het verleden samenhangen met de huidige cijfers. In de eerste twee kwartalen van 2020 wordt bij een viertal variabelen een sterke afwijking gemeten in de internetreeks. We zien het zorggebruik sterk afnemen. Het gaat hier om de variabelen huisartscontact, tandartsbezoek en specialistcontact. Bij de variabele ervaren gezondheid zien we juist een sterke toename. Voor deze variabelen geldt de aanname van de samenhang met het verleden niet meer en passen we het model aan. Dit is gesignaleerd doordat de gestandaardiseerde innovaties waardes aannemen die (absoluut) veel groter zijn dan 2 (zie vorige paragraaf).

Het tijdreeksmodel leent informatie uit het verleden via zowel de trend \(L_t\) als het seizoen \(S_t\). De verspreiding van corona kan zowel invloed hebben op de trend als het seizoenspatroon. Aangezien het niet mogelijk is het effect op het seizoenspatroon apart te schatten nemen we in het model aan dat corona alleen effect heeft op de dynamiek van de trend. In het model is de trend gemodelleerd als een smooth-trendmodel (paragraaf 2.1), gegeven door 
\begin{eqnarray*}
L_t & = & L_{t-1} + R_{t-1} \\
R_t & = & R_{t-1} + \eta_t^{\rm R}, \mbox{ waarbij}
\end{eqnarray*}
\begin{align*}
& \eta_t^{\rm R} \sim {\cal N}(0,\sigma_{\rm R}^2)  \\
& {\rm Cov}(\eta_t^{\rm R},\eta_{t'}^{\rm R}) = 0, \mbox{voor $t\neq t'$.}
\end{align*}

De trend is opgebouwd uit een niveau \(L_t\) en een hellingsparameter \(R_t\) (de slope), waarbij de variantie \(\sigma_{\rm R}^2\) van de storingsterm \(\eta_t^{\rm R}\) van de slope constant is over de tijd. De variantie bepaald de flexibiliteit van de trend. Voor de variabelen over het zorggebruik en de ervaren gezondheid maken we de trend in het model flexibeler door de variantie \(\sigma_{\rm R}^2\) met een tijdsafhankelijke factor \(f_t \geq 1\) te vermenigvuldigen:
\begin{align*}
& \eta_t^{\rm R} \sim {\cal N}(0, f_t \sigma_{\rm R}^2)  \\
& {\rm Cov}(\eta_t^{\rm R},\eta_{t'}^{\rm R}) = 0, \mbox{voor $t\neq t'$.}
\end{align*}

De kwartalen waar \(f_t >1\) gekozen wordt en de waarden van \(f_t\) bij die kwartalen worden bepaald aan de hand van de gestandaardiseerde innovaties. Hierbij zorgen we ervoor dat de gestandaardiseerde innovaties in alle kwartalen niet te ver buiten het interval (-2,2) vallen, zodat de normaliteit van de innovaties gegarandeerd blijft. Ook proberen we de variantie van de slope storingstermen \(f_t \sigma_{\rm R}^2 \) zo min mogelijk aan te passen (dat wil zeggen, waarden voor \(f_t\) zo klein mogelijk houden), zodat het model nog informatie kan lenen uit het verleden. Merk op dat het aanpassen van \(\sigma_{\rm R}^2\) in kwartaal \(t\) invloed heeft op de storingsterm in kwartaal \(t\) en daarom op de slope pas in kwartaal \(t+1\) en op de trend pas in kwartaal \(t+2\). Er is dus een vertraging van twee kwartalen in het effect op de uitkomsten na aanpassing van \(\sigma_{\rm R}^2\). 

Na analyse van de gestandaardiseerde innovaties bij de variabelen over het zorggebruik de variantie \(\sigma_{\rm R}^2\) aangepast van kwartaal 3 in 2019 t/m kwartaal 2 in 2020 en bij ervaren gezondheid vanaf kwartaal 2 in 2019 t/m kwartaal 2 in 2020. Bij de andere variabelen is het niet nodig om de slope flexibeler te maken, omdat de gestandaardiseerde innovaties daar geen aanleiding toe geven. Ook in latere kwartalen van 2020 blijkt het niet nodig te zijn om de slope flexibeler te maken. Tabel 1 geeft een overzicht van de modelaanpassingen per doelvariabele. 

Na deze aanpassingen worden er bij de meeste doelvariabelen geen modelveronderstellingen verworpen. Bij het berekenen van de eerste twee kwartalen wordt alleen bij huisartscontact en dagelijks roken de veronderstelde normaliteit verworpen. Bij deze variabelen blijken in enkele kwartalen nog uitbijters voor te komen, maar de gestandaardiseerde innovaties in deze kwartalen vallen net buiten het interval (-2,2). Bij huisartscontact gaat het om kwartaal 3 in 2016 en kwartaal 1 in 2020. Bij dagelijks roken gaat het om kwartaal 2 in 2019. Bij het berekenen van de kwartalen 3 en 4 wordt er geen enkele modelveronderstelling meer verworpen.