---
title: "Uge 3 - løsningsforslag"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Øvelse 1

## Hent pakker

```{r}
library("tidyverse")
library("broom")
```



## Hent data

```{r}
fev <- read.csv(url("http://biostatistics.dk/puff/data/fev.csv"), header=TRUE)
head(fev)
```

## Exploratory data analysis 

```{r}
boxplot(fev$FEV ~ fev$Smoke)
```



## Model (Spgm A + B + C)

```{r}
fev$Smoke <- factor(fev$Smoke)
fev$Smoke <- relevel(fev$Smoke, "Yes")
model1 <- lm(FEV ~ Smoke, data=fev)
model1 %>% tidy()
model1 %>% confint()

# 
fev$tal <- 1:6
fev$tal <- factor(fev$tal)

lm(FEV ~ tal, data=fev) %>% tidy()
```

Jeg er 95% sikker på, at intervallet [0.49 ; 0.93] indeholder den sande forskel i middelværdier 
for FEV mellem børn, der ryger og børn der ikke ryger


```{r}
model2 <- lm(FEV ~ Smoke + Age, data=fev)
model2 %>% tidy()
model2 %>% confint()
```


Jeg er 95% sikker på, at intervallet [0.05 ; 0.36] indeholder den sande forskel i middelværdi
for FEV mellem børn, der ryger og børn der ikke ryger, når der er korrigeret/justeret for at børnene har forskellige aldre.


## Prædiktioner (spgm D + E)

12-årig, der ikke ryger

```{r}
0.158 + 0.209 + 12*0.231
```


12-årig, der ryger

```{r}
0.158 + 0.209*0 + 12*0.231
```

```{r}
DF <- data.frame(Smoke=c("No", "Yes"), Age=c(12, 12))
predict(model2, newdata=DF, interval="prediction")
```


## Vekselvirkningsanalyse (spgm F + G)

```{r}
model3 <- lm(FEV ~ Smoke*Age, data=fev)
model3 %>% tidy()
model3 %>% confint()

```


Modelkontrol

```{r}
library("MESS")
model4 <- lm(FEV ~ Smoke*Age, data=fev)
wallyplot(model4)
residual_plot(model4)

```

## Prædiktion (spgm H + I)

```{r}
DF <- data.frame(Smoke=c("No", "Yes"), Age=c(12, 12))
predict(model3, newdata=DF, interval="prediction")
```

Vekselvirkningen lader ikke til at bidrage til modellen, så ville afrapportere resultaterne fra den simplere model uden vekselvirkning.

# AMH


Den primære hypotese er at se på afhængigheden af alder på hormonet AMH, samt effekten af, hvorvidt kvinden indtager P-piller.

Hent data

```{r}
load(url("http://biostatistics.dk/puff/data/amh.rda"))
```

Er interesseret i at modellere amh som funktion af de øvrige forklarende variable.


Opstiller først en model, der indeholder de forklarende variable.

```{r}
m1 <- lm(amh ~ alder + bmi + rygning + ppiller, data=amh)
```

Laver dernæst modelkontrol for at checke om modellen lader til at holde

```{r}
residualplot(m1)
```

Det kunne godt se ud som variansheterogenitet. Lad os forsøge med at logaritmetransformere udfaldet for at se om det bliver bedre.

```{r}
m2 <- lm(log(amh) ~ alder + bmi + rygning + ppiller, data=amh)
residualplot(m2)
```

Det ser lidt bedre ud, så lad os fortsætte med denne model. Bemærk, at resultaterne af de efterfølgende udregninger alle er på logaritmeskala.

# Estimater

```{r}
m2 %>% tidy()
confint(m2)
```

De primære formål var at undersøge alder og p-pillers indflydelse på AMH. Af analysen kan vi se, at både p-piller og alder har negative koefficient, altså vil amh falde med alderen og hvis man bruger p-piller. Effekten af både alder og p-piller er ret kraftig, og eksempelvis lader brugen af p-piller til at have ca. samme effekt som en forskel på 4-5 år.

Vi kan desuden undersøge om der kunne være en vekselvirkning mellem p-piller og alder: er det sådan, at effekten af p-piller på hormonet er afhænging af kvindens alder?

```{r}
m3 <- lm(log(amh) ~ alder*ppiller + bmi + rygning , data=amh)
m3 %>% tidy()
m3 %>% confint()
```

Her ses det ud fra sidste linje, at de parameteren svarende til vekselvirkning er ret tæt på 0, og den derfor kan udelades. 95% konfidensintervallet går fra -0.05 til 0.003.

Vi kan også se, at vekselvirkningen kan droppes ved at bruge `drop()`.

```{r}
drop1(m3, test="F")
```

Igen ses det, at vekselvirkningen ikke bidrager en masse til modelleringen.

## Tilbagetransformation

Alle udregningerne ovenfor er baseret på de logaritmerede hormonværdier. Derfor vil vi gerne tilbagetransformere resultaterne til den oprindelige skala. Estimaterne bliver

```{r}
exp(coef(m2))
```

Og da disse er tilbagetransformeret fra logaritmen skal de fortolkes som *forholdsmæssige* ændringer. Værdien 0.923 for alder betyder, at for hvert år ekstra skal vi *gange* en faktor 0.923 på hormonniveauet. Det svarer til at hormonniveauet hvert år falder med 1-0.923 = 0.077 = 7.7%.

For p-piller har vi et estimat på 0.700. Det betyder, at kvinderne der bruger p-piller har et niveau af hormonet, der er 30% under niveauet for kvinder, der ikke bruger p-piller.

Vi kan desuden udregne et 95% konfidensinterval for parametrene på den tilbagetransformerede skala:

```{r}
exp(confint(m2))
```

Vi er altså eksempelvis 95% sikre på, at intervallet [ 0.61 ; 0.80 ] indeholder det sande *forhold* i niveau af hormonet amh mellem gruppen af kvinder, der bruger p-piller og gruppen, der ikke gør.




