---
title: "Vitalkapacitet"
output: html_document
---

```{r message=FALSE}
library(mosaic)
```

I denne øvelse analyseres et datasæt vedrørende vitalkapacitet, som er den maksimale luftmængde, der kan udåndes efter en maksimal indånding.
Datasættet kan hentes via dette [link](http://asta.math.aau.dk/dan/static/datasets?file=vitcap.txt), og gemmes
i samme mappe som denne Rmd-fil, så du kan køre tingene helt på din egen 
computer. Alternativt kan du også bare indlæse data direkte fra web-adressen,
hvor det så hentes hver gang du kører kommandoen (f.eks. ved at trykke "knit"):
```{r}
vitcap <- read.delim("http://asta.math.aau.dk/dan/static/datasets?file=vitcap.txt")
vitcap <- vitcap %>% mutate(exposure = as.factor(exposure))
```


I datasættet findes variablen `vital.capacity`, som er blevet målt for 84 arbejdere i cadmiumindustrien.

En yderligere variabel er faktoren `exposure`, som har 3 niveauer, der angiver, hvor længe den enkelte arbejder har været eksponeret for cadmium:

- A: Ingen eksponering
- B: Mindre end 10 år
- C: Mere end 10 år

Derudover inderholder datasættet også variablen `age`, der angiver den enkelte arbejders alder.
(Og yderligere dummy-variablene `z1` og `z2`, som er beskrevet i ekstraopgaven, hvor de benyttes.)

## Opgave 1 (ANOVA-analyse af exposures effekt på vitalkapacitet)

Lav en analyse, der undersøger effekten af faktoren `exposure` på responsen `vital.capacity`.
Dvs. en lineær model på formen:
$$
y = \alpha + \beta_1 z_1 + \beta_2 z_2 + \varepsilon
$$
hvor $z_1$ og $z_2$ er dummy variable der indikerer om man er hhv. i gruppe B og C.

Du bør som minimum:

- Lave et relevant boxplot.
```{r}
# gf_boxplot(...)
```

- Lave summary af modellen.
```{r}
# model1 <- lm(..., data = vitcap)
# summary(model1)
```

- Opskriv prædiktionsligningen
$$
\hat{y} = \hat\alpha + \hat\beta_1 z_1 + \hat\beta_2 z_2
$$
hvor du indsætter tallene fra det tidligere summary for parameterestimaterne.

- Hvad er den forventede forskel i lungekapacitet mellem personer i gruppe A og C?
- Gennemgå F-testen for ingen effekt af `exposure` (findes i slutningen af summary):
    + Hvad er nulhypotesen?
    + Hvad er test-statistikken?
    + Hvad er p-værdien?
    + Hvilken konklusion træffer vi om den samlede effekt af exposure på vitalkapaciteten?


## Opgave 2 (Additiv model for alders og exposures effekt på vitalkapacitet)

Vi udvider nu analysen til også at omfatte arbejdernes alder (`age`) som en forklarende variabel for responsen `vital.capacity`.

Tilpas en model hvor `age` og `exposure` indgår additivt:
$$
y = \alpha + \beta_1 z_1 + \beta_2 z_2 + \beta_3 x + \varepsilon
$$
hvor $x$ er alder og $z_1$ og $z_2$ angiver hhv. gruppe B og C som før.
```{r}
# model2 <- lm(...)
# summary(model2)
```

Gennemgå detaljerne for denne model:

- Illustrer modellen grafisk med f.eks. `plotModel()`. (Tip: Vær opmærksom på, at `plotModel` forventer, at den definerede model har den kontinuerte variabel `age` til at stå først i formeludtrykket. Så hvis man vil plotte regressionslinjerne for den additive model, så skal man bruge formalutrykket `vital.capacity ~ age + exposure`, når man tilpasser sin model, og IKKE `vital.capacity ~ exposure + age`. De to formeludtryk repræsenterer ganske vist præcis den samme additive model, men den bagvedliggende kode i funktionen `plotModel` virker desværre kun efter hensigten, hvis `age` står først.)

- Med udgangspunkt i outputtet af et `summmary` for modellen:
    + Indsæt talværdier i den samlede prædiktionsligning.
$$
\hat{y} = \hat\alpha + \hat\beta_1 z_1 + \hat\beta_2 z_2 + \hat\beta_3 x
$$
    + Opskriv prædiktionsligningerne for hver af de tre grupper (hint: indsæt værdierne for dummyvariable og lad led med 0 udgå).
$$
A: \quad \hat{y} = ...
$$
$$
B: \quad \hat{y} = ...
$$
$$
C: \quad \hat{y} = ...
$$
    + Hvad er den forventede vitalkapacitet for en arbejder på 40 år i gruppe B?
    + Hvor meget forventes vitalkapaciteten at ændre sig årligt for en arbejder
      i gruppe C? Hvad med for dem i gruppe A og B?
    + Er effekten af alder signifikant? Gennemgå herunder teststørrelse og p-værdi og hvordan disse hænger sammen.
    + Er der en signifikant forskel mellem gruppe A og B? Hvad med A og C?

## Opgave 3 (Interaktionsmodel for alders og exposures effekt på vitalkapacitet)

Vi udvider nu analysen til også at tillade en interaktion mellem alder og exposures effekt på vitalkapaciteten.

Tilpas en model hvor `age` og `exposure` indgår med interaktionseffekt
$$
y = \alpha + \beta_1 z_1 + \beta_2 z_2 + \beta_3 x + \beta_4 z_1 x + \beta_5 z_2 x + \varepsilon
$$
hvor $x$ er alder og $z_1$ og $z_2$ angiver hhv. gruppe B og C som før.
```{r}
# model3 <- lm(...)
# summary(model3)
```

Gennemgå detaljerne for denne model:

- Illustrer modellen grafisk med f.eks. `plotModel()` som før.
- Med udgangspunkt i outputtet af et `summmary` for modellen:
    + Indsæt talværdier i den samlede prædiktionsligning.
$$
\hat{y} = \hat\alpha + \hat\beta_1 z_1 + \hat\beta_2 z_2 + \hat\beta_3 x + \hat\beta_4 z_1 x + \hat\beta_5 z_2 x
$$
    + Opskriv prædiktionsligningerne for hver af de tre grupper (hint: indsæt værdierne for dummyvariable og lad led med 0 udgå).
$$
A: \quad \hat{y} = ...
$$
$$
B: \quad \hat{y} = ...
$$
$$
C: \quad \hat{y} = ...
$$
    + Hvad er den forventede vitalkapacitet for en arbejder på 40 år i gruppe B?
    + Hvor meget forventes vitalkapaciteten at ændre sig årligt for en arbejder
      i gruppe C? Hvad med for dem i gruppe A og B?
    + Er der signifikant forskel i effekten af alder mellem gruppe A og B? Hvad med mellem A og C?
    + Brug `anova()` til at undersøge om der samlet set er en signifikant interaktion mellem effekten af exposure og alder på vitalkapaciteten.

## Opfølgende spørgsmål:

Hvordan kan det være, at det i opgavens første model (uden `age` som forklarende variabel) ser ud som om, at vi kan tillade os at slå de tre `exposure`-grupper sammen, mens det modellerne med `age` som forklarende variabel ikke ser sådan ud?


## Ekstraopgave (hvis du har tid)
Datasættet indeholder også dummy-variabler for faktoren `exposure`:

- `z1=1` hvis `exposure=B` ellers 0
- `z2=1` hvis `exposure=C` ellers 0

Betragt følgende to modeller:
```{r}
model4 <- lm(vital.capacity ~ age*z2, data = vitcap)
model5 <- lm(vital.capacity ~ age*z1 + age*z2, data = vitcap)
```

- Brug en F-test til at vise, at der ikke er nogen signifikant forskel mellem `model4` og `model5` (Hint: bruge `anova()`).
- Giv en fortolkning af forskellen mellem de to modeller.

I forbindelse med det sidste spørgsmål kan du lave følgende plot:

```{r}
gf_point(vital.capacity ~ age, data = vitcap, color = ~ factor(z2)) %>% 
  gf_smooth(method = "lm", se = F) %>% 
  gf_labs(color = "Gruppe C")
```