6. Enkelvoudige lineaire regressie

1 Breast cancer dataset

• subset van studie https://doi.org/10.1093/jnci/djj052

• 32 borstkanker patiënten met een estrogen recepter positieve tumor die tamoxifen chemotherapy behandeling ondergaan. Variabelen:

  • grade: histologische graad van tumor (graad 1 vs 3),
  • node: status van de lymfe knopen (0: niet aangetast, 1: aantasting en verwijdering van de lymfe knopen),
  • size: grootte van tumor in cm,
  • ESR1 en S100A8 gen expressie in tumor biopsy (via microarray technologie)

brca <- read_csv("https://raw.githubusercontent.com/statOmics/sbc20/master/data/breastcancer.csv")
brca

# A tibble: 32 x 10
   sample_name filename     treatment    er grade  node  size   age  ESR1 S100A8
   <chr>       <chr>        <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
 1 OXFT_209    gsm65344.ce… tamoxifen     1     3     1   2.5    66 1939.  207. 
 2 OXFT_1769   gsm65345.ce… tamoxifen     1     1     1   3.5    86 2752.   37.0
 3 OXFT_2093   gsm65347.ce… tamoxifen     1     1     1   2.2    74  379. 2364. 
 4 OXFT_1770   gsm65348.ce… tamoxifen     1     1     1   1.7    69 2532.   23.6
 5 OXFT_1342   gsm65350.ce… tamoxifen     1     3     0   2.5    62  141. 3219. 
 6 OXFT_2338   gsm65352.ce… tamoxifen     1     3     1   1.4    63 1495.  108. 
 7 OXFT_2341   gsm65353.ce… tamoxifen     1     1     1   3.3    76 3406.   14.0
 8 OXFT_1902   gsm65354.ce… tamoxifen     1     3     0   2.4    61 2813.   68.4
 9 OXFT_1982   gsm65355.ce… tamoxifen     1     1     0   1.7    62  950.   74.2
10 OXFT_5210   gsm65356.ce… tamoxifen     1     3     0   3.5    65 1053.  182. 
# … with 22 more rows

• Om didactische redenen verwijderen we eerst 3 outliers in de S100A8 expressie data.
• In deze studie kan dit echter niet worden verantwoord
• Later in de lessen laten we zien hoe je goed met alle data omgaat.

brca %>%
  ggplot(aes(x="", y=S100A8)) +
  geom_boxplot() +
  xlab("") +
  ylab("S100A8 expressie")

---

library(GGally)
brcaSubset <- brca %>%
  filter(S100A8<2000)
brcaSubset[,-(1:4)] %>% ggpairs()

1.1 Associatie tussen ESR1 en S100A8 expressie

• ESR1 in \(\pm\) 75% van borstkankertumoren.

  • Expressie van ER-gen positief voor behandeling: tumor reageert op hormoontherapie
  • Tamoxifen interageert met ER en moduleert genexpressie.

• Eiwitten van de S100-familie zijn vaak gedisreguleerd bij kanker

• S100A8 expressie onderdrukt immuunsysteem in tumor en creëert inflamatoir milieu die kankergroei promoot.

• Interesse in associatie tussen ESR1 en S100A8 expressie.

1. pipe dataset naar ggplot
2. selecteer data ggplot(aes(x=ESR1,y=S100A8))
3. voeg punten toe met geom_point()
4. voeg een “smooth line” toe geom_smooth()

brcaSubset %>%
  ggplot(aes(x=ESR1,y=S100A8)) +
  geom_point() +
  geom_smooth()

2 Lineaire Regressie

• Statistische methode om relatie tussen 2 reeksen observaties \((X_i, Y_i)\), bekomen voor onafhankelijke subjecten \(i = 1, ..., n\), te beschrijven.

• Gen expressie voorbeeld

  • Response Y : S100A8 expressie
  • Predictor X: ESR1 expressie

1. pipe dataset naar ggplot
2. selecteer data ggplot(aes(x=ESR1,y=S100A8))
3. voeg punten toe met geom_point()
4. voeg een “smooth line” toe geom_smooth()
5. voeg een rechte toe geom_smooth() met method = "lm" (linear model). (We zetten se = FALSE om geen puntgewijze betrouwbaarheidsintervallen weer te geven)

brcaSubset %>%
  ggplot(aes(x = ESR1,y = S100A8)) +
  geom_point() +
  geom_smooth(se = FALSE, col = "grey") +
  geom_smooth(method = "lm", se = FALSE)

2.1 Model

• Voor vaste \(X\), heeft \(Y\) niet noodzakelijke dezelfde waarde

\[\text{observation = signal + noise}\]

\[Y_i=g(X_i)+\epsilon_i\]

• We definiëren \(g(x)\) als het verwachte resultaat voor subjects met \(X_i=x\)

\[E[Y_i|X_i=x]=g(x)\]

Daarom is \(\epsilon_i\) gemiddeld 0 voor subjects met dezelfde \(X_i\): \[E[\epsilon_i|X_i]=0\]

2.2 Lineaire regressie

• Om accurate en interpreteerbare resultaten te bekomen, kiest men \(g(x)\) vaak als een lineaire functie met ongekende parameters.

\[E(Y|X=x)=\beta_0 + \beta_1 x\]

onbekend intercept \(\beta_0\) en helling \(\beta_1\).

• Lineair model legt een assumptie op de verdeling van \(X\) en \(Y\), die incorrect kan zijn.

• Efficiënte data-analyse: benut alle observaties om iets te leren over verwachte uitkomst bij \(X=x\).

2.3 Gebruik

• Predictie: wanneer \(Y\) ongekend is, maar \(X\) wel, kunnen we \(Y\) voorspellen op basis van \(X\) \[E(Y|X=x)=\beta_0 + \beta_1 x\]

• Associatie: biologische relatie tussen variabele \(X\) en continue meting \(Y\) beschrijven.

• Intercept: \(E(Y|X=0)=\beta_0\)

• Slope: \[\begin{eqnarray*} E(Y|X=x+\delta)-E(Y|X=x)&=&\beta_0 + \beta_1 (x+\delta) -\beta_0-\beta_1 x\\ &=& \beta_1\delta \end{eqnarray*}\]

\(\beta_1:\) verschil in gemiddelde uitkomst voor subjecten die verschillen in één eenheid van de predictor \(X\).

3 Parameterschatting

• • Kleinste kwadraten techniek (Least squares)

brcaSubset %>%
  ggplot(aes(x = ESR1, y = S100A8)) +
  geom_point() +
  geom_smooth(se = FALSE, col = "grey") +
  geom_smooth(method = "lm", se = FALSE)

• Parameters \(\beta_0\) en \(\beta_1\) zijn ongekend

• Parameters schatten op basis van beperkte steekproef

• Best passende lijn

  • Punt op regressielijn voor een gegeven \(x_i\): \((x_i, \beta_0 + \beta_1 x_i)\) zo dicht mogelijk bij \((x_i, y_i)\)
  • Kies \(\beta_0\) en \(\beta_1\) zodat de som tussen voorspelde en waargenomen punten zo klein mogelijk wordt.

\[SSE=\sum_{i=1}^n (y_i-\beta_0-\beta_1 x_i)^2=\sum_{i=1}^n e_i^2\]

• Met \(e_i\) de residuen: de verticale afstanden van de observaties tot de gefitte regressierechte

3.1 Schatters die SSE minimaliseren

\[\hat{\beta_1}= \frac{\sum\limits_{i=1}^n (y_i-\bar y)(x_i-\bar x)}{\sum\limits_{i=1}^n (x_i-\bar x_i)^2}=\frac{\mbox{cor}(x,y)s_y}{s_x} \]

\[\hat{\beta_0}=\bar y - \hat{\beta}_1 \bar x \]

• Merk op dat de helling van de kleinste kwadratenlijn evenredig is met de correlatie tussen de uitkomst en de verklarende variabele.

Geschatte lineaire regressiemodel laat toe om:

• verwachte uitkomst te voorspellen voor subjecten met een gegeven waarde \(x\) voor de predictor: \[\text{E} [ Y | X = x]=\hat{\beta}_0+\hat{\beta}_1x\]

• na te gaan hoeveel uitkomst gemiddeld verschilt tussen 2 groepen subjecten met een verschil van \(\delta\) eenheden in de verklarende variabele:

\[\text{E}\left[Y|X=x+\delta\right]-\text{E}\left[Y|X=x\right]= \hat{\beta}_1\delta\]

3.1.1 Borstkanker voorbeeld

lm1 <- lm(S100A8 ~ ESR1, brcaSubset)
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

\[E(Y|X=x)=208.47-0.059 x\]

• De verwachte S100A8-expressie is gemiddeld 59 eenheden lager voor patiënten met een ESR1-expressieniveau die 1000 eenheden hoger ligt.

• Verwachte S100A8 expressieniveau voor patiënten met een ESR1 expressieniveau van 2000:
  \[208.47-0.059\times 2000=89.94\]

• Verwachte S100A8 expressieniveau voor patiënten met een ESR1 expressieniveau van 4000:
  \[208.47-0.059\times 4000=-28.58\]

• Let op voor extrapolatie! (Veronderstelling van lineariteit kan men enkel nagaan binnen het bereik van de data).

4 Statistische besluitvorming (statistische inferentie)

Om besluiten te kunnen trekken over lineaire regressiemodel \[E(Y|X)=\beta_0+\beta_1 X\]

moeten we weten:

• Hoe de least squares parameter schatters variëren van steekproef tot steekproef, en

• En hoe is dit onder de nulhypothese dat er geen associatie is tussen predictor en response

• Noodzaak aan statistisch model!

• Modelleer de verdeling van \(Y\) gegeven \(X\) expliciet: \(f_{Y|X}(y)\)

4.1 Modelleer verdeling van Y?

1. Naast lineariteit hebben we nog aannames nodig!
2. Onafhankelijkheid: waarnemingen \((X_1,Y_1), ..., (X_n,Y_n)\) zijn gemaakt voor n onafhankelijke subjecten (Is vereist om de variantie te schatten)
3. Homoscedasticiteit of gelijke varianties: waarnemingen variëren met gelijk gemiddelde rond de regressielijn.
   • Residuen \(\epsilon_i\) hebben gelijke variantie voor elke \(X_i=x\)
   • \(\text{var}(Y\vert X=x) = \sigma^2\) voor elke \(X=x\)
   • \(\sigma\) wordt de residuele standaarddeviatie genoemd.
4. Normaliteit: de residuen \(\epsilon_i\) zijn normaal verdeeld

• Uit 2, 3 en 4 volgt dat \[\epsilon_i \text{ i.i.d.} N(0,\sigma^2).\]

• Als we ook steunen op eerste veronderstelling van lineariteit: \[Y_i\vert X_i\sim N(\beta_0+\beta_1 X_i,\sigma^2),\]

• 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)\]

4.2 Hoge spreiding op \(X\) verbetert de precisie

\[\sigma^2_{\hat{\beta}_1}=\frac{\sigma^2}{\sum\limits_{i=1}^n (X_i-\bar X)^2}\]

• Conditionele variantie (\(\sigma^2\)) is niet gekend
• Schatten d.m.v. gemiddelde van die kwadratische afwijkingen rond de regressierechte
• mean squared error (MSE)

\[\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). - deel door \(n-2\)

Na schatting van \(\sigma^2\) bekomen we volgende standaard errors:

\[\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}}\]

• Opnieuw toetsen en betrouwbaarheidsintervallen o.b.v. \[T=\frac{\hat{\beta}_k-\beta_k}{SE(\hat{\beta}_k)} \text{ with } k=1,2.\]

• Als aan alle aannames is voldaan volgt \(T\) een t-verdeling met n-2 vrijheidsgraden.

• Als geen normaliteit maar wel onafhankelijk, lineariteit en homoscedasticiteit en grote dataset \[\rightarrow \text{Centrale Limietstelling}\]

4.2.1 Borstkanker voorbeeld

• Negatieve associatie tussen S100A8 en ESR1 gen expressie.

• Veralgemeen effect in steekproef naar populatie met behulp van het betrouwbaarheidsinterval op 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

• Negatief verband is significant op het 5% significantieniveau.

4.3 Hypothese test

• Vertaal de onderzoeksvraag “Is er een associatie tussen de S100A8- en ESR1-genexpressie?” naar de parameters in het model.

• Onder nulhypothese geen associatie tussen expressie van beide genen: \[H_0: \beta_1=0\]

• Onder alternatieve hypothese is er een associatie tussen beide genen: \[H_1: \beta_1\neq0\]

• Test statistiek \[T=\frac{\hat{\beta}_1-0}{SE(\hat{\beta}_k)}\]

• Onder \(H_0\) volgt de statistiek een t-verdeling met n-2 vrijheidsgraden.

4.3.1 Brca dataset

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

• Associatie tussen de S100A8 en ESR1 genexpressie extreem significant is (p<<0.001).
• Maar eerst moeten we alle assumpties controleren!
• Anders kunnen de conclusies o.b.v. de statistische test en het BI onjuist zijn.

5 Nagaan van modelveronderstellingen

• Onafhankelijkheid: design
• Lineariteit: besluitvorming geen zin als model niet lineair is
• Homoscedasticiteit: besluitvorming/p-waarde is niet betrouwbaar als de data niet homoscedastisch zijn
• Normaliteit: besluitvorming/p-waarde is niet betrouwbaar als de data niet normaal verdeeld zijn in kleine steekproeven

5.1 Lineariteit

brcaSubset %>%
  ggplot(aes(x = ESR1, y = S100A8)) +
  geom_point() +
  geom_smooth(se = FALSE, col = "grey") +
  geom_smooth(method = "lm", se = FALSE)

5.2 Residuplots

• Afwijkingen van lineariteit echter makkelijker opgespoord d.m.v. een residuplot. (Zeker als er later meer variabelen zijn in het lineaire model)

• verklarende variabele of predicties \(\hat\beta_0+\hat\beta_1 x\) op de \(X\)-as

• de residuen op de \(Y\)-as \[e_i=y_i-\hat{g}(x_i)=y_i-\hat\beta_0-\hat\beta_1\times x_i,\]

plot(lm1)

5.3 Homoscedasticiteit (gelijkheid van variantie)

• Residuen en kwadratische residu’s dragen informatie over residuele variabiliteit.

• Associatie met de verklarende variabelen \(\rightarrow\) indicatie van heteroscedasticiteit.

• Scatterplot van of \(e_i\) versus \(x_i\) of predicties \(\hat \beta_0+ \hat \beta_1 x_i\).

• Scatterplot van gestandardiseerd residu versus \(x_i\) of predicties.

5.4 Normaliteit

• Indien voldoende gegevens, schatters normaal verdeeld zelfs wanneer observaties niet Normaal verdeeld zijn: centrale limiet stelling

• Wat ‘voldoende observaties’ zijn, hangt af van hoe goed de verdeling op de Normale lijkt.

• Aanname is dat uitkomsten Normaal verdeeld zijn bij vaste waarden van de verklarende variabelen. \[Y_i\vert X_i\sim N(\beta_0+\beta_1X_i,\sigma^2)\]

• QQ-plot van response Y is heel misleidend.

• QQ-plot van de residuen \(e_i\)

plot(lm1, which=2)

6 Afwijkingen van Modelveronderstellingen

• Transformatie van onafhankelijke veranderlijke wijzigt de verdeling van Y bij gegeven X niet:

  • helpt niet om normaliteit of homoscedasticiteit te bekomen
  • helpt wel om lineariteit te bekomen wanneer er normaliteit en homoscedasticiteit is
  • Vaak ook hogere orde termen: \(X^2\), \(X^3\), … \[Y_i=\beta_0+\beta_1X_i+\beta_2X_i^2+ ... + \epsilon_i\]

• Transformatie van response Y kan helpen om normaliteit en homoscedasticiteit te bekomen.

  • \(\sqrt(Y)\), \(\log(Y)\), 1/Y, …

6.1 Borstkanker voorbeeld

Problemen met

• heteroscedasticiteit

• mogelijkse afwijking van normaliteit (scheefheid naar rechts)

• negatieve concentratievoorspellingen die theoretisch niet mogelijk zijn

• niet-lineairiteit

• treedt veelal op bij concentratie en intensiteitsmetingen

• Deze zijn vaak log-normaal verdeeld (normale verdeling na log-transformatie)

• In Figuur 6.3 eveneens een soort exponentiële trend

• In de genexpressie literatuur wordt veelal gebruik gemaak van \(\log_2\) transformatie

• gen-expressie op log-schaal proportionele verschillen op de originele schaal

brca %>%
  ggplot(aes(x = ESR1, y = S100A8)) +
  geom_point() +
  geom_smooth()

brca %>%
  ggplot(aes(
    x = ESR1 %>% log2,
    y = S100A8 %>% log2)
    ) +
  geom_point() +
  geom_smooth()

lm2<-lm(S100A8 %>% log2 ~ ESR1 %>% log2, brca)
plot(lm2)

summary(lm2)

Call:
lm(formula = S100A8 %>% log2 ~ ESR1 %>% log2, data = brca)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.94279 -0.66537  0.08124  0.68468  1.92714 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     23.401      1.603   14.60 3.57e-15 ***
ESR1 %>% log2   -1.615      0.150  -10.76 8.07e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.026 on 30 degrees of freedom
Multiple R-squared:  0.7942,    Adjusted R-squared:  0.7874 
F-statistic: 115.8 on 1 and 30 DF,  p-value: 8.07e-12

confint(lm2)

                  2.5 %    97.5 %
(Intercept)   20.128645 26.674023
ESR1 %>% log2 -1.921047 -1.308185

6.1.1 Interpretatie 1

Een patiënt met een ESR1 expressie die 1 eenheid op de \(\log_2\) schaal hoger ligt dan dat van een andere patiënt heeft gemiddeld gezien een expressie-niveau van het S100A8 gen dat 1.61 eenheden lager ligt (95% BI [-1.92,-1.31]).

Crossectionele studie: enkel uitspraken over verschillen tussen patiënten!

\[\log_2 \hat\mu_1=23.401 -1.615 \times \text{logESR}_1,\text{ } \log_2 \hat\mu_2=23.401 -1.615 \times \text{logESR}_2 \] \[\log_2 \hat\mu_2-\log_2 \hat\mu_1= -1.615 (\log_2 \text{ESR}_2-\log_2 \text{ESR}_1) = -1.615 \times 1 = -1.615\]

6.1.2 Interpretatie 2

Model op log-schaal: bij terugtransformatie verkrijgen we het geometrische gemiddelde

\[\begin{eqnarray*} \sum\limits_{i=1}^n \frac{\log x_i}{n}&=&\frac{\log x_1 + \ldots + \log x_n}{n}\\\\ &\stackrel{(1)}{=}&\frac{\log(x_1 \times \ldots \times x_n)}{n}=\frac{\log\left(\prod\limits_{i=1}^n x_i\right)}{n}\\\\ &\stackrel{(2)}{=}&\log \left(\sqrt[\leftroot{-1}\uproot{2}\scriptstyle n]{\prod\limits_{i=1}^n x_i}\right) \end{eqnarray*}\]

• Populatiegemiddelden \(\mu\) dus geschat a.d.h.v. geometrisch gemiddelden.
• Logaritmische transformatie is een monotoon: we kunnen betrouwbaarheidsintervallen berekend op log-schaal terugtransformeren!

2^lm2$coef[2]

ESR1 %>% log2 
    0.3265519 

2^-lm2$coef[2]

ESR1 %>% log2 
       3.0623 

2^-confint(lm2)[2,]

   2.5 %   97.5 % 
3.786977 2.476298 

Een patiënt met een ESR1 expressie die 2 keer zo hoog is als die van een andere patiënt, zal gemiddeld een S100A8-expressie hebben die 3.06 keer lager is (95% BI [2.48,3.79]).

\[\log_2 \hat\mu_1=23.401 -1.615 \times \text{logESR}_1,\text{ } \log_2 \hat\mu_2=23.401 -1.615 \times \text{logESR}_2 \] \[\log_2 \hat\mu_2-\log_2 \hat\mu_1= -1.615 (\log_2 \text{ESR}_2-\log_2 \text{ESR}_1) \] \[\log_2 \left[\frac{\hat\mu_2}{\hat\mu_1}\right]= -1.615 \log_2\left[\frac{ \text{ESR}_2}{\text{ESR}_1}\right] \] \[\frac{\hat\mu_2}{\hat\mu_1}=\left[\frac{ \text{ESR}_2}{\text{ESR}_1}\right]^{-1.615}=2^ {-1.615} =0.326\] or \[\frac{\hat\mu_1}{\hat\mu_2}=2^{1.615} =3.06\]

6.1.3 Interpratie 3

Een patiënt met een ESR1 expressie die 1% hoger is dan die van een andere patiënt zal gemiddeld een expressieniveau voor het S100A8 gen hebben dat ongeveer -1.61% lager is (95% BI [-1.92,-1.31])%.

\[\log_2 \hat\mu_1=23.401 -1.615 \times \text{logESR}_1,\text{ } \log_2 \hat\mu_2=23.401 -1.615 \times \text{logESR}_2 \] \[\log_2 \hat\mu_2-\hat\log_2 \mu_1= -1.615 (\log_2 \text{ESR}_2-\log_2 \text{ESR}_1) \] \[\log_2 \left[\frac{\hat\mu_2}{\hat\mu_1}\right]= -1.615 \log_2\left[\frac{ \text{ESR}_2}{\text{ESR}_1}\right] \] \[\frac{\hat\mu_2}{\hat\mu_1}=\left[\frac{ \text{ESR}_2}{\text{ESR}_1}\right]^{-1.615}=1.01^ {-1.615} =0.984 \approx -1.6\%\]

Dit geldt voor lage tot matige waarden van \(\beta_1\): \[-10<\beta_1<10 \rightarrow 1.01^{\beta_1} -1 \approx \frac{\beta_1}{100}.\]

7 Besluitvorming over gemiddelde uitkomst

• Regressie model kan ook worden gebruikt voor predictie

• Besluitvorming te doen over de gemiddelde uitkomst bij een gegeven waarde \(x\), m.a.w. \[\hat{g}(x)= \hat{\beta}_0 + \hat{\beta}_1 x\]

• \(\hat{g}(x)\) een schatter van het conditionele gemiddelde \(E[Y\vert X=x]\)

• Parameterschatters Normale verdeeld en onvertekend \(\rightarrow\) schatter \(\hat{g}(x)\) ook Normaal verdeeld en onvertekend.

\[\text{SE}_{\hat{g}(x)}=\sqrt{MSE\left\{\frac{1}{n}+\frac{(x-\bar X)^2}{\sum\limits_{i=1}^n (X_i-\bar X)^2}\right\}}.\]

\[T=\frac{\hat{g}(x)-g(x)}{SE_{\hat{g}(x)}}\sim t_{n-2}\]

• Gemiddelde uitkomst en betrouwbaarheidsintervallen op de gemiddelde uitkomst in R via de predict(.) functie.
• newdata argument: predictorwaarden (x-waarden) voor het berekenen van gemiddelde uitkomsten
• interval="confidence" argument om betrouwbaarheidsintervallen te bekomen.
• Zonder newdata argument wordt de gemiddelde uitkomsten berekend voor alle predictorwaarden van de dataset.

grid <- 140:4000
g <- predict(
  lm2,
  newdata = data.frame(ESR1 = grid),
  interval = "confidence")
head(g)

       fit      lwr      upr
1 11.89028 10.76082 13.01974
2 11.87370 10.74721 13.00019
3 11.85724 10.73370 12.98078
4 11.84089 10.72028 12.96151
5 11.82466 10.70696 12.94237
6 11.80854 10.69372 12.92336

Merk op dat we de nieuwe data die we gespecificeerd hebben voor de ESR1 expressie niet moeten transformeren omdat we het model fitten met de lmfunctie en de transformatie hebben gespecificeerd binnen die functie met behulp van het pipe commando!

brca %>% ggplot(
  aes(
    x = ESR1 %>% log2,
    y = S100A8 %>% log2)
  ) +
  geom_point() +
  geom_smooth(method = "lm")

7.1 Terugtransformatie

newdata<-data.frame(cbind(grid,2^g))
brca %>%
  ggplot(aes(x = ESR1, y = S100A8)) +
  geom_point() +
  geom_line(aes(x=grid,y=fit),newdata) +
  geom_line(aes(x=grid,y=lwr),newdata,color="grey") +
  geom_line(aes(x=grid,y=upr),newdata,color="grey")

8 Predictie-intervallen

-We kunnen ook een voorspelling doen voor de locatie van een nieuwe waarneming die zou worden verzameld in een nieuw experiment voor een patiënt met een bepaalde waarde voor hun ESR1-expressie

• Het is belangrijk op te merken dat dit experiment nog moet worden uitgevoerd. We willen dus de niet-waargenomen individuele expressiewaarde voor een nieuwe patiënt voorspellen.

• Voor een nieuwe onafhankelijke waarneming \(Y^*\) \[ Y^* = g(x) + \epsilon^* \] met \(\epsilon^*\sim N(0,\sigma^2)\) en \(\epsilon^*\) onafhankelijk van de waarnemingen in de steekproef \(Y_1,\ldots, Y_n\).

• We voorspellen een nieuwe log-S100A8 voor een patiënt met een gekend log2-ESR1-expressieniveau x \[ \hat{y}(x)=\hat{\beta}_0+\hat{\beta}_1 \times x \]

• De geschatte gemiddelde uitkomst en voorspelling voor een nieuwe waarneming zijn gelijk.

• Maar hun steekproef verdelingen zijn anders!

  • Onzekerheid over de geschatte gemiddelde uitkomst \(\leftarrow\) onzekerheid over de geschatte modelparameters \(\hat\beta_0\) en \(\hat\beta_1\).
  • Onzekerheid over nieuwe waarneming \(\leftarrow\) onzekerheid over geschat gemiddelde en extra onzekerheid omdat de nieuwe waarneming zal afwijken rond het gemiddelde!

\[\text{SE}_{\hat{Y}(x)}=\sqrt{\hat\sigma^2+\hat\sigma^2_{\hat{g}(x)}}=\sqrt{MSE\left\{1+\frac{1}{n}+\frac{(x-\bar X)^2}{\sum\limits_{i=1}^n (X_i-\bar X)^2}\right\}}.\]

\[\frac{\hat{Y}(x)-Y}{\text{SE}_{\hat{Y}(x)}}\sim t_{n-2}\]

• Merk op dat een predictie-interval (PI) een verbeterde versie is van een referentie-interval wanneer de modelparameters onbekend zijn: onzekerheid over modelparameters + t-verdeling.

p <- predict(
  lm2,
  newdata = data.frame(ESR1 = grid),
  interval="prediction")

head(p)

       fit      lwr      upr
1 11.89028 9.510524 14.27004
2 11.87370 9.495354 14.25205
3 11.85724 9.480288 14.23419
4 11.84089 9.465324 14.21646
5 11.82466 9.450461 14.19886
6 11.80854 9.435698 14.18138

preddata<-data.frame(
  grid = grid%>%log2,
  p)
brca %>% ggplot(aes(x=ESR1%>%log2,y=S100A8%>%log2)) +
  geom_point() +
  geom_smooth(method="lm") +
     geom_line(aes(x=grid,y=lwr),preddata,color="blue") +
  geom_line(aes(x=grid,y=upr),preddata,color="blue")

preddata<-data.frame(cbind(grid,2^p))
brca %>% ggplot(aes(x = ESR1, y = S100A8)) +
  geom_point() +
  geom_line(
    aes(x = grid,y = fit),
    newdata) +
  geom_line(
    aes(x = grid, y = lwr),
    newdata,
    color = "grey") +
  geom_line(
    aes(x = grid, y = upr),
    newdata,
    color = "grey") +
  geom_line(
    aes(x = grid, y = lwr),
    preddata,
    color = "blue") +
  geom_line(
    aes(x = grid,y = upr),
    preddata,
    color = "blue")

8.1 NHANES voorbeeld

• Vergelijk referentie-interval voor cholesterolgehalte met predictie interval.

• Referentie-interval

library(NHANES)
fem <- NHANES %>%
  filter(Gender=="female"&!is.na(DirectChol))

2^(
  fem %>%
  pull(DirectChol) %>%
  log2 %>%
  mean +
    c(-1,1) *
    qnorm(0.975) *
    (fem %>%
      pull(DirectChol) %>%
      log2 %>%
      sd)
  )

[1] 0.8361311 2.4397130

• Predictie interval

lmChol <- lm(DirectChol %>% log2 ~ 1, data=fem)
predInt <- predict(
  lmChol,
  interval="prediction",
  newdata=data.frame(noPred=1)
  )
round(2^predInt,2)

   fit  lwr  upr
1 1.43 0.84 2.44

Merk op dat het voorspellingsinterval bijna gelijk is aan het referentie-interval voor de grote steekproef. We konden de parameters inderdaad heel precies schatten.

We zullen hetzelfde doen voor een kleine steekproef van 10 patiënten.

• Referentie interval

set.seed(1)
fem10 <- NHANES %>%
  filter(Gender=="female"&!is.na(DirectChol)) %>%
  sample_n(size=10)

2^(
  fem10 %>%
    pull(DirectChol) %>%
    log2 %>%
    mean +
      c(-1,1) *
      qnorm(0.975) *
      (fem10 %>%
        pull(DirectChol) %>%
        log2 %>%
        sd)
  )

[1] 0.8976012 2.2571645

Het referentie-interval is veel smaller dan in de grote steekproef.

• Predictie interval

lmChol10 <- lm(DirectChol %>% log2 ~ 1, data = fem10)
predInt10 <- predict(
  lmChol10,
  interval = "prediction",
  newdata = data.frame(noPred=1)
  )
round(2^predInt10, 2)

   fit  lwr  upr
1 1.42 0.81 2.49

• Merk op dat het PI nu onzekerheid meeneemt in parameterschatters (gemiddelde en standaard error). En dat het interval veel breder wordt! Dit is hier vooral belangrijk voor de bovengrens omdat we de gegevens terug hebben getransformeerd!

• Het interval is bijna net zo breed als dat gebaseerd op de grote steekproef.

• Bij kleine steekproeven is het erg belangrijk om met deze extra onzekerheid rekening te houden.

9 Kwadratensommen en Anova-tabel

9.1 Totale kwadratensom

\[\text{SSTot} = \sum_{i=1}^n (Y_i-\bar{Y})^2.\]

• SStot kan worden gebruikt om de variantie van de marginale verdeling van de respons te schatten.

• In dit hoofdstuk hebben we ons gefocused op de conditionele verdeling \(f(Y\vert X=x)\).

• We weten dat MSE een goede schatting is van de variantie van de conditionele verdeling van \(Y\vert X=x\).

9.2 Kwadratensom van de regressie SSR

\[\text{SSR} = \sum_{i=1}^n (\hat{Y}_i - \bar{Y})^2 = \sum_{i=1}^n (\hat{g}(x_i) - \bar{Y})^2.\]

• Maat voor de afwijking tussen de predicties op de geschatte regressierechte en het steekproefgemiddelde van de uitkomsten.

• Een andere interpretatie: verschil tussen twee modellen

  • Geschatte model \(\hat{g}(x)=\hat\beta_0+\hat\beta_1x\)
  • Geschatte model zonder predictor (enkel intercept): \(g(x)=\beta_0\) \(\rightarrow\) \(\beta_0\) zal gelijk zijn aan \(\bar{Y}\).

• SSR meet de grootte van het effect van de predictor

9.3 Kwadratensom van de fouten

\[ \text{SSE} = \sum_{i=1}^n (Y_i-\hat{Y}_i )^2 = \sum_{i=1}^n \left\{Y_i-\hat{g}\left(x_i\right)\right\}^2.\]

• Hoe kleiner de SSE, hoe beter het model fit.

• Kleinste kwadraten techniek!

---

We kunnen aantonen dat SST kan worden ontbonden in \[\begin{eqnarray*} \text{SSTot} &=& \sum_{i=1}^n (Y_i-\bar{Y})^2 \\ &=& \sum_{i=1}^n (Y_i-\hat{Y}_i+\hat{Y}_i-\bar{Y})^2 \\ &=& \sum_{i=1}^n (Y_i-\hat{Y}_i)^2+\sum_{i=1}^n(\hat{Y}_i-\bar{Y})^2 \\ &=& \text{SSE }+\text{SSR} \end{eqnarray*}\]

• De totale variabiliteit in de gegevens (SSTot) wordt gedeeltelijk verklaard door het regressieverband (SSR).
• Variabiliteit die we niet kunnen verklaren met het regressiemodel is de residuele variabiliteit (SSE).

9.4 Determinatie-coëfficiënt

\[ R^2 = 1-\frac{\text{SSE}}{\text{SSTot}}=\frac{\text{SSR}}{\text{SSTot}}.\]

• Fractie van de totale variabiliteit in de steekproef-uitkomsten die verklaard wordt door het geschatte regressieverband.

• Grote \(R^2\) indicatie dat model potentieel tot goede predicties kan leiden (kleine SSE)

• Slechts in beperkte mate indicatief voor de p-waarde van de test \(H_0:\beta_1=0\) vs \(H_1:\beta_1\neq0\).

  • p-waarde sterk beïnvloed door SSE en steekproefgrootte \(n\), maar niet door SSTot
  • De determinatiecoëfficiënt \(R^2\) wordt door SSE en SSTot bepaald, maar niet door de steekproefgrootte n.

• Model met lage \(R^2\) blijft wel nuttig om associatie te bestuderen, zolang het de associatie correct modelleert!

9.4.1 Borstkanker voorbeeld

summary(lm2)

Call:
lm(formula = S100A8 %>% log2 ~ ESR1 %>% log2, data = brca)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.94279 -0.66537  0.08124  0.68468  1.92714 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     23.401      1.603   14.60 3.57e-15 ***
ESR1 %>% log2   -1.615      0.150  -10.76 8.07e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.026 on 30 degrees of freedom
Multiple R-squared:  0.7942,    Adjusted R-squared:  0.7874 
F-statistic: 115.8 on 1 and 30 DF,  p-value: 8.07e-12

9.5 F-Testen in het enkelvoudig lineair regressiemodel

• Kwadratensommen zijn basis voor \(F\)-tests \[ F = \frac{\text{MSR}}{\text{MSE}}\]

met \(\text{MSR} = \frac{\text{SSR}}{1} \text{ and } \text{MSE} = \frac{\text{SSE}}{n-2}.\)

• MSR wordt de gemiddelde kwadratensom van de regressie genoemd.

• noemers 1 en \(n-2\) zijn de vrijheidsgraden van SSR en SSE.

• onder \(H_0: \beta_1=0\) volgt de teststatistiek \[H_0:F = \frac{\text{MSR}}{\text{MSE}} \sim F_{1,n-2},\]

• F-test is altijd twee-zijdig! \(H_1:\beta_1\neq 0\) \[ p = P_0\left[F\geq f\right]=1-F_F(f;1,n-2)\]

summary(lm2)

Call:
lm(formula = S100A8 %>% log2 ~ ESR1 %>% log2, data = brca)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.94279 -0.66537  0.08124  0.68468  1.92714 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     23.401      1.603   14.60 3.57e-15 ***
ESR1 %>% log2   -1.615      0.150  -10.76 8.07e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.026 on 30 degrees of freedom
Multiple R-squared:  0.7942,    Adjusted R-squared:  0.7874 
F-statistic: 115.8 on 1 and 30 DF,  p-value: 8.07e-12

9.6 Anova Tabel

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Regressie	vrijheidsgraden SSR	SSR	MSR	f-statistiek	p-waarde
Error	vrijheidsgraden SSE	SSE	MSE

anova(lm2)

Analysis of Variance Table

Response: S100A8 %>% log2
              Df  Sum Sq Mean Sq F value   Pr(>F)    
ESR1 %>% log2  1 121.814 121.814   115.8 8.07e-12 ***
Residuals     30  31.559   1.052                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

10 Dummy variabelen

• Lineaire regressiemodel voor het vergelijken van twee gemiddelden.

• Borstkanker: verschil is in gemiddelde leeftijd van de patiënten met onaangetaste lymfeknopen en patiënten waarvan lymfeknopen werden verwijderd.

• Hiervoor definiëren we eerst een \(dummy\) variabele \[x_i = \left\{ \begin{array}{ll} 1 & \text{aangetaste lymfeknopen} \\ 0 & \text{onaangetaste lymfeknopen} \end{array}\right.\]

• groep met \(x_i=0\) wordt de referentiegroep genoemd.

• Het regressiemodel blijft ongewijzigd, \[Y_i = \beta_0 + \beta_1 x_i +\epsilon_i\] met \(\epsilon_i \text{ iid } N(0,\sigma^2)\)

Gezien \(x_i\) slechts twee waarden kan aannemen, is het eenvoudig om het regressiemodel voor beide waarden van \(x_i\) afzonderlijk te bekijken: \[ \begin{array}{lcll} Y_i &=& \beta_0 +\epsilon_i &\text{onaangetaste lymfeknopen} (x_i=0) \\ Y_i &=& \beta_0 + \beta_1 +\epsilon_i &\text{ aangetaste lymfeknopen} (x_i=1) . \end{array}\]

Dus \[\begin{eqnarray*} E\left[Y_i\mid x_i=0\right] &=& \beta_0 \\ E\left[Y_i\mid x_i=1\right] &=& \beta_0 + \beta_1, \end{eqnarray*}\]

waaruit direct de interpretatie van \(\beta_1\) volgt: \[ \beta_1 = E\left[Y_i\mid x_i=1\right]-E\left[Y_i\mid x_i=0\right]\]

\(\beta_1\) is dus het gemiddelde verschil in leeftijd tussen patiënten met aangetaste lymfeknopen en patiënten met onaangetaste lymfeknopen (referentiegroep).

Met de notatie \(\mu_0= E\left[Y_i\mid x_i=0\right]\) en \(\mu_1= E\left[Y_i\mid x_i=1\right]\) wordt dit \[\beta_1 = \mu_1-\mu_0.\]

Er kan aangetoond worden dat \[\begin{array}{ccll} \hat\beta_0 &=& \bar{Y}_1&\text{ (steekproefgemiddelde in referentiegroep)} \\ \hat\beta_1 &=& \bar{Y}_2-\bar{Y}_1&\text{(schatter van effectgrootte)} \\ \text{MSE} &=& S_p^2 . \end{array}\]

De testen voor \(H_0:\beta_1=0\) vs. \(H_1:\beta_1\neq0\) kunnen gebruikt worden voor het testen van de nulhypothese van de two-sample \(t\)-test, \(H_0:\mu_1=\mu_2\) t.o.v. \(H_1:\mu_1\neq\mu_2\).

brca$node <- as.factor(brca$node)
t.test(age~node,brca,var.equal=TRUE)

    Two Sample t-test

data:  age by node
t = -2.7988, df = 30, p-value = 0.008879
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -15.791307  -2.467802
sample estimates:
mean in group 0 mean in group 1 
       59.94737        69.07692 

lm3 <- lm(age~node, brca)
summary(lm3)

Call:
lm(formula = age ~ node, data = brca)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.9474  -5.3269   0.0526   5.3026  18.0526 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   59.947      2.079  28.834  < 2e-16 ***
node1          9.130      3.262   2.799  0.00888 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.063 on 30 degrees of freedom
Multiple R-squared:  0.207, Adjusted R-squared:  0.1806 
F-statistic: 7.833 on 1 and 30 DF,  p-value: 0.008879

plot(lm3)

brca %>%
  ggplot(
    aes(
      x = node %>%
        as.factor,
      y = age)
      ) +
  geom_boxplot()

1. Simulate 9 datasets with the same number of observations as the brca dataset from a normal distribution with the same standard deviation as in the original data. Store the data of all simulations in a data frame
2. Plot the simulated data using the ggplot function
3. Add a boxplot layer
4. Use facet_wrap to make a separate plot for simulated dataset
5. Change label of y-axis

set.seed(354)
sim_df <- data.frame(
  node = rep(brca$node, 9),
  iid = rnorm(9 * nrow(brca), sd = sigma(lm3)),
  sim = rep(1:9, each = 32)
  )
sim_df %>%
  ggplot(aes(x = node, y=iid)) +
  geom_boxplot() +
  facet_wrap(.~sim) +
  ylab(paste0("iid N(0,",round(sigma(lm3)^2,2),")"))

10.1 Observationele study

• We kunnen echter niet besluiten dat oudere personen een hoger risico hebben op aantasting van de lymfeknopen ten gevolge van hun leeftijd.

• Mogelijks confounding: geen randomisatie \(\rightarrow\) groepen patiënten met aangetaste lymfeknopen en niet-aangetaste lymfeknopen kunnen nog in andere karateristieken van elkaar verschillen.

• Enkel besluiten dat er een associatie is tussen de lymfeknoop status en de leeftijd.

• Het is dus niet noodzakelijkerwijs een causaal verband!

---

• Is ook zo voor lineair model voor de \(\log_2\)-S100A8-expressie.

• Aangezien we de ESR1-expressie niet experimenteel vast hebben kunnen leggen, kunnen we niet besluiten dat een hogere ESR1-expressie de S100A8-expressie doet verlagen.

• Enkel besluiten dat er een negatieve associatie is.

• Om impact van gen te bestuderen op andere genen: knockout mutanten generenen in het labo

---

Home