Als de gegevens representatief zijn voor de populatie kan men in de regressiecontext eveneens aantonen dat de kleinste kwadraten schatters voor het intercept en de helling onvertekend zijn, m.a.w
\[E[\hat \beta_0]=\beta_0 \text{ en } E[\hat \beta_1]=\beta_1\]Het feit dat de schatters gemiddeld (over een groot aantal vergelijkbare studies) niet afwijken van de waarden in de populatie, impliceert niet dat ze niet rond die waarde variëren. Om inzicht te krijgen hoe dicht we de parameterschatters bij het werkelijke intercept \(\beta_0\) en de werkelijke helling \(\beta_1\) mogen verwachten, wensen we bijgevolg ook haar variabiliteit te kennen.
In de borstkanker dataset hebben we een negatieve associatie geobserveerd tussen de S100A8 en ESR1 gen expressie. Net zoals in Hoofdstuk 5 is het op basis van de puntschatters voor de helling niet duidelijk of dat verband werkelijk voorkomt in de populatie of indien we het verband door toeval hebben geobserveerd in de dataset. De schatting van de helling is immers onnauwkeurig en zal variëren van steekproef tot steekproef. Het resultaat van een data-analyse is dus niet interpreteerbaar zonder die variabiliteit in kaart te brengen.
Om de resultaten uit de steekproef te kunnen veralgemenen naar de populatie zullen we in deze context eveneens inzicht nodig hebben op de verdeling van de parameterschatters. Om te kunnen voorspellen hoe de parameterschatters variëren van steekproef tot steekproef enkel en alleen op basis van slechts één steekproef zullen we naast de onderstelling van
bijkomende aannames moeten maken over de verdeling van de gegevens, met name
Uit 2, 3 en 4 volgt dus dat de residuen \(\epsilon_i\) onafhankelijk zijn en dat ze allen eenzelfde Normale verdeling volgen
\[\epsilon_i \sim N(0,\sigma^2).\]Als we ook steunen op de veronderstelling van lineariteit weten we dat de originele observaties conditioneel op \(X\) eveneens Normaal verdeeld zijn
\[Y_i\sim N(\beta_0+\beta_1 X_i,\sigma^2),\]met een gemiddelde dat varieert in functie van de waarde van de onafhankelijke variabele \(X_i\).
Verder kan men aantonen dat onder deze aannames
\[\sigma^2_{\hat{\beta}_0}=\frac{\sum\limits_{i=1}^n X^2_i}{\sum\limits_{i=1}^n (X_i-\bar X)^2} \times\frac{\sigma^2}{n} \text{ en } \sigma^2_{\hat{\beta}_1}=\frac{\sigma^2}{\sum\limits_{i=1}^n (X_i-\bar X)^2}\]en dat de parameterschatters eveneens normaal verdeeld zijn
\[\hat\beta_0 \sim N\left(\beta_0,\sigma^2_{\hat \beta_0}\right) \text{ en } \hat\beta_1 \sim N\left(\beta_1,\sigma^2_{\hat \beta_1}\right)\]Merk op dat de onzekerheid op de helling af zal nemen wanneer er meer observaties zijn en/of wanneer de observaties meer gespreid zijn. Voor het opzetten van een experiment kan dit belangrijke informatie zijn. Uiteraard wordt de precisie ook beïnvloed door de grootte van de variabiliteit van de observaties rond de rechte, \(\sigma^2\), maar dat heeft een onderzoeker meestal niet in de hand.
De conditionele variantie (\(\sigma^2\)) is echter niet gekend en is noodzakelijk voor de berekening van de variantie op de parameterschatters. We kunnen \(\sigma^2\) echter ook schatten op basis van de observaties. Zoals beschreven in Hoofdstuk 4 kunnen we de variatie van de uitkomsten rond hun conditionele gemiddelde beschrijven d.m.v. de afwijkingen tussen de observaties \(y_i\) en hun (geschatte) gemiddelde \(\hat{g}(x)=\hat{\beta}_0+\hat{\beta}_1x_i\), de residu’s. Het gemiddelde van die residu’s is echter altijd 0 omdat positieve en negatieve residu’s mekaar opheffen. Bijgevolg levert het gemiddelde residu geen goede maat op voor de variatie en is het beter om naar kwadratische afwijkingen \(e_i^2\) te kijken. Net zoals de steekproefvariantie een goede schatter was voor de variantie (Sectie 4.3.2), zal in de regressiecontext het gemiddelde van die kwadratische afwijkingen rond de regressierechte opnieuw een goede schatter zijn voor \(\sigma^2\). Deze schatter wordt in de literatuur ook wel de mean squared error (MSE) genoemd.
\[\hat\sigma^2=MSE=\frac{\sum\limits_{i=1}^n \left(y_i-\hat\beta_0-\hat\beta_1\times x_i\right)^2}{n-2}=\frac{\sum\limits_{i=1}^n e^2_i}{n-2}.\]Voor het bekomen van deze schatter steunen we op onafhankelijkheid (aanname 2) en homoscedasticiteit (aanname 3). Merk op dat we bij deze schatter niet delen door het aantal observaties \(n\), maar door \(n-2\). Hierbij corrigeren we voor het feit dat voor de berekening van MSE 2 vrijheidsgraden worden gespendeerd aan het schatten van het intercept en de helling.
Na het schatten van MSE kunnen we \(\sigma^2\) door MSE vervangen zodat schatters worden bekomen voor de variantie en standard error op de schatters van model parameters,
\[\text{SE}_{\hat{\beta}_0}=\hat\sigma_{\hat{\beta}_0}=\sqrt{\frac{\sum\limits_{i=1}^n X^2_i}{\sum\limits_{i=1}^n (X_i-\bar X)^2} \times\frac{\text{MSE}}{n}} \text{ en } \text{SE}_{\hat{\beta}_1}=\hat\sigma_{\hat{\beta}_1}=\sqrt{\frac{\text{MSE}}{\sum\limits_{i=1}^n (X_i-\bar X)^2}}\]Analoog als in Hoofdstuk 5 kunnen we opnieuw toetsen en betrouwbaarheidsintervallen construeren op basis van de teststatistieken
\[T=\frac{\hat{\beta}_k-\beta_k}{SE(\hat{\beta}_k)} \text{ met } k=0,1.\]Als aan alle aannames is voldaan dan volgen deze statistieken \(T\) een t-verdeling met n-2 vrijheidsgraden. Wanneer niet is voldaan aan de veronderstelling van normaliteit maar wel aan lineariteit, onafhankelijkheid en homoscedasticiteit dan kunnen we voor inferentie opnieuw beroep doen op de centrale limietstelling die zegt dat de statistiek T bij benadering een standaard Normaal verdeling zal volgen wanneer het aantal observaties voldoende groot is.
In de borstkanker dataset hebben we een negatieve associatie geobserveerd tussen de S100A8 en ESR1 gen expressie. We kunnen het effect in de steekproef nu veralgemenen naar de populatie toe door een betrouwbaarheidsinterval te bouwen voor de helling:
\[[\hat\beta_1 - t_{n-2,\alpha/2} \text{SE}_{\hat\beta_1},\hat\beta_1 + t_{n-2,\alpha/2} \text{SE}_{\hat\beta_1}]\]confint(lm1)
## 2.5 % 97.5 %
## (Intercept) 149.84639096 267.09649989
## ESR1 -0.08412397 -0.03440378
Op basis van de R-output bekomen we een 95% betrouwbaarheidsinterval voor de helling \(-0.084,-0.034\). Gezien nul niet in het interval ligt weten we eveneens dat de negatieve associatie statistisch significant is op het 5% significantieniveau.
Anderzijds kunnen we ook een formele hypothesetoets uitvoeren. Onder de nulhypothese veronderstellen we dat er geen associatie is tussen de expressie van beide genen:
\[H_0: \beta_1=0\]en onder de alternatieve hypothese is er een associatie tussen beide genen:
\[H_1: \beta_1\neq0\]Met de test statistiek
\[T=\frac{\hat{\beta}_1-0}{SE(\hat{\beta}_k)}\]kunnen we de nulhypothese falsifiëren. Onder \(H_0\) volgt de statistiek een t-verdeling met n-2 vrijheidsgraden.
Deze tweezijdige test is geïmplementeerd in de standaard output van R.
summary(lm1)
##
## Call:
## lm(formula = S100A8 ~ ESR1, data = brcaSubset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -95.43 -34.81 -6.79 34.23 145.21
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 208.47145 28.57207 7.296 7.56e-08 ***
## ESR1 -0.05926 0.01212 -4.891 4.08e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 59.91 on 27 degrees of freedom
## Multiple R-squared: 0.4698, Adjusted R-squared: 0.4502
## F-statistic: 23.93 on 1 and 27 DF, p-value: 4.078e-05
De test geeft weer dat de associatie tussen de S100A8 en ESR1 genexpressie extreem significant is (p<<0.001). Als de nulhypothese waar is en als aan alle voorwaarden is voldaan dan is er een kans van 4 op 100000 om een helling te vinden die minstens even extreem is door toeval. Het is bijgevolg heel onwaarschijnlijk om dergelijke associatie te observeren in een steekproef wanneer de nulhypothese waar is.
Vooraleer we een conclusie trekken is het echter belangrijk dat we alle aannames verifiëren omdat de statistische test en de betrouwbaarheidsintervallen anders incorrect zijn.