---
title: "Practicum 2: Stokstaartjes"
output:
  html_document:
    code_download: yes
    highlight: tango
    number_sections: yes
    theme: cosmo
    toc: yes
    toc_float: yes
  pdf_document:
    toc: yes
  word_document:
    toc: yes
---

Stokstaartjes, een mangoestensoort die voorkomt in zuidelijk Afrika, leven in groepen. Elke groep kent zijn dominant vrouwtje, die grotendeels het monopolie heeft op de voortplanting binnen de groep: ze is de moeder van meer dan tachtig procent van de jongen die hun eerste maand overleven. Hoe ze dit monopolie bewerktstelligt, is echter minder goed begrepen. Een mogelijke verklaring ligt in het feit dat het dominante vrouwtje ongeveer een maand voor ze bevalt de andere vrouwtjes terroriseert, wat voor hen vaak uitmondt in een tijdelijke verbanning uit de groep. In een onderzoek tracht men na te gaan of deze verbanning gerelateerd is met factoren die een nadelige invloed hebben op de voortplanting. 
De studiepopulatie waarop het onderzoek betrekking heeft, bedraagt individuen uit 15 willekeurig gekozen groepen stokstaartjes. In een eerste luik van het experiment werd het gewichtsverschil bekeken van 15 stokstaartjes die verbannen werden uit de groep. In een tweede luik van het experiment werden op gelijkaardige tijdstippen als het vertrek en de terugkomst van de verbannen vrouwtjes het gewichtsverschil gemeten van 15 andere stokstaartjes uit dezelfde regio, maar die in hun eigen groep konden blijven. Deze laatste groep van individuen worden aangeduid in de dataset met `verbannen == 0` .


#Gewichtsverschil met en zonder verbanning
Test op het 5% significantieniveau de onderstelling dat de gemiddelde gewichtsverandering bij vrouwtjes die verbannen waren uit de groep **verschilt van** vrouwtjes die in de groep konden blijven. Noteer tevens duidelijk uw conclusie.

# R libraries inlezen

```{r}
library(ggplot2)
library(dplyr)
library(tidyr)
library(readr)
```

# Vraag 1 Dataset inlezen

De dataset `stokstaartjesDat.txt` bestaat uit $60$ observaties van de volgende variabelen:

 - `gewicht` het gewicht van het stokstaartje
 - `tijdstip` tijdstip van de meting. `"voor"` indien gewogen werd voor de verbanning, anders `"na"`
 - `verbannen` $1$ indien verbannen, $0$ indien niet.

```{r}
#In het tekstbestand zijn er drie kolomnamen gegeven op de eerste lijn; gewicht, tijdstip en verbannen. Er zijn echter vier waarden vanaf de tweede lijn, waarbij de eerste waarde een rij nummer is. Dit zal dus voor problemen zorgen bij het inlezen. Om dit probleem te verhelpen gebruiken we skip = 1 om de eerste rij met de kolomnamen over te slaan. Daarna geven we zelf kolomnamen door met col_names. Met col_types vermelden we specifiek welk type variabele iedere kolom is. Met col_select selecteren we dan de kolommen die we wel degelijk nodig hebben.

stokstaartjes <- read_delim("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum3/stokstaartjesDat.txt", 
                            delim = "\ " , 
                            skip = 1, 
                            col_names = c("rij_nr","gewicht","tijdstip","verbannen"), 
                            col_types = list(col_integer(), col_double(), col_factor(), col_factor()), 
                            col_select = c("gewicht","tijdstip","verbannen"))

stokstaartjes <- stokstaartjes %>% 
  mutate(
	  index = factor(c(rep(1:15,2),rep(16:30,2)))
  )
```

# Vraag 2 Data exploratie
We willen visueel nagaan of er verschil zichtbaar is in het gewicht voor en na verbanning van de groep. Maak een ggplot die dit toont voor de groep die verbannen is en voor de groep die niet verbannen is. 

Deze plot moet dus de individuele metingen van het gewicht tonen per tijdstip en de paring van de data tonen aan de hand van de index variabele. De plot kan je splitsen per groep (variabele verbannen) aan de hand van de facet_wrap layer (zie ?facet_wrap). Sla deze figuur op in `lineplot`.

```{r}
lineplot <- ... %>% ggplot(...) + 
  ... + # om lineplot te maken
  facet_wrap(...) + 
  theme_bw()
lineplot
```

# Vraag 3 Data exploratie
Zet de dataset om in wide format en maak een kolom `verschil` aan dat het verschil bevat tussen de na en voor metingen. Sla dit dataframe op in `stokstaartjes_wide`.

Maak een boxplot `boxplot_verschil` van het verschil in gewicht van de stokstaartjes die in de groep bleven en van de stokstaartjes die verbannen zijn. Geef de boxplots van de groepen verschillende kleuren. Is er een verschil tussen de groepen?


```{r}
stokstaartjes_wide <- stokstaartjes %>% 
  pivot_wider(id_cols = c("verbannen","index"), #id_cols zijn de kolommen die iedere observatie uniek identificeren
              names_from = ..., 
              values_from = ...) %>% 
  mutate(...)

boxplot_verschil <- stokstaartjes_wide %>% ggplot()+ 
  ... + #boxplot layer met argument outlier.shape = NA
  ... + #observaties toevoegen, gebruik een jitter
  theme_bw() 
boxplot_verschil
```

# Vraag 4 Test (Meerkeuzevraag)
Welke test zullen we uitvoeren om deze onderzoeksvraag op te lossen? 

Ter herinnering, de onderzoeksvraag luidt als volgt: Test op het 5% significantieniveau de onderstelling dat de gemiddelde gewichtsverandering bij vrouwtjes die verbannen waren uit de groep verschilt van vrouwtjes die in de groep konden blijven

1. ongepaarde, tweezijdige t-test 
2. gepaarde, tweezijdige t-test
3. ongepaarde, eenzijdige t-test
4. gepaarde, eenzijdige t-test

# Vraag 5 Test 

We voeren een ongepaarde t-test (two-sample t-test) uit omdat de groepen die vergeleken worden onafhankelijk zijn van elkaar. De groepen zijn onafhankelijk omdat de vergelijking wordt gemaakt met een willekeurig ander vrouwtje uit de regio dat in haar groep kon blijven. 

De onderzoekers zijn geïnteresseerd in een verschil in gewichtsverandering, ongeacht de richting. Het is bijvoorbeeld ook mogelijk dat de stokstaartjes in groep het eten moeten delen, en dus mogelijks nog meer gewicht verliezen dan de verbannen individuen. Daarom zullen we tweezijdig testen.


# Vraag 6 Test: Voorwaarden
We veronderstellen opnieuw onafhankelijkheid van de gegevens op basis van het studiedesign. 
We zullen de assumptie van normaliteit binnen elke groep nagaan aan de hand van een QQ-plot, gebruik facet_wrap (variabele verbannen) om de plot te splitsen per groep. Sla deze plot op in `qq_plot_verschil`.

```{r}
qq_plot_verschil <- stokstaartjes_wide %>% ggplot(...) +
  ... +
  ... +
  facet_wrap(...) + 
  theme_bw()

qq_plot_verschil
```

# Vraag 7 Test: Voorwaarden
De voorwaarde van normaliteit lijkt vrij goed voldaan te zijn. De punten liggen dicht rond de rechte bij de verbannen stokstaartjes. Bij de stokstaartjes die bleven is de afwijking al iets groter. We kunnen opnieuw eens vergelijken met gesimuleerde data onder de nulhypothese. De afwijkingen die we in onze qqplot zien lijken niet uitzonderlijk te zijn.

```{r}
set.seed(1406)
nobs <- sum(stokstaartjes_wide$verbannen == 0)
n_simulations <- 9

data.frame(
  y = c(rnorm(nobs*n_simulations,
              mean = stokstaartjes_wide$verschil[stokstaartjes_wide$verbannen == 0],
              sd = sd(stokstaartjes_wide$verschil[stokstaartjes_wide$verbannen == 0])
             )
      ),
  label = rep(paste0("sim",1:n_simulations)),
              each = nobs) %>%
  ggplot(aes(sample = y)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~ label, scales = "fixed")+ 
  ylab("Sample Quantiles") +
  xlab("Theoretical Quantiles") +
  ggtitle("Simulatie QQ Plots Verbannen Groep") + 
  theme_bw()


data.frame(
  y = c(rnorm(nobs*n_simulations,
              mean = stokstaartjes_wide$verschil[stokstaartjes_wide$verbannen == 1],
              sd = sd(stokstaartjes_wide$verschil[stokstaartjes_wide$verbannen == 1])
             )
      ),
  label = rep(paste0("sim",1:n_simulations)),
              each = nobs) %>%
  ggplot(aes(sample = y)) +
  geom_qq() +
  geom_qq_line() +
  facet_wrap(~ label, scales = "fixed")+ 
  ylab("Sample Quantiles") +
  xlab("Theoretical Quantiles") +
  ggtitle("Simulatie QQ Plots Niet Verbannen Groep") + 
  theme_bw()
```

# Vraag 8 Test: Voorwaarden
Om mogelijks power te winnen, zullen we nagaan of de varianties gelijk zijn tussen beide groepen. In het geval dat ze voldoende gelijk zijn, kunnen we er immers voor kiezen om de t-test met gepoolde steekproefvariantie uit te voeren in plaats van de t-test met Welch-benadering. Bekijk de eerder gemaakte boxplot `boxplot_verschil` nogmaals om de voorwaarden na te gaan.

```{r}
boxplot_verschil
```

# Vraag 9 Test: Voorwaarden
Het is vaak niet eenvoudig om te beoordelen of de varianties sterk van mekaar verschillen of niet. Om een beter idee te krijgen, kunnen we eens een aantal boxplots simuleren met dezelfde steekproefgrootte en in de veronderstelling dat de varianties gelijk zijn. 
Daarnaast bestaat er ook een test in R die test of de varianties al dan niet gelijk zijn. De nulhypothese van deze test is dat de varianties gelijk zijn, de alternatieve hypothese is dat ze verschillend zijn.

```{r}
set.seed(52)
par(mfrow=c(3,3), mar=c(3,2,1,1))
sd1<-sd(stokstaartjes_wide$verschil) 
mean1<-mean(stokstaartjes_wide %>% filter(verbannen == 1) %>% pull(verschil))
mean2<-mean(stokstaartjes_wide %>% filter(verbannen == 0) %>% pull(verschil))
for(i in 1:9){
  x<-rnorm(15, mean=mean1, sd=sd1)
  y<-rnorm(15, mean=mean2, sd=sd1)
  boxplot(x,y)
}

var.test(stokstaartjes_wide$verschil~stokstaartjes_wide$verbannen)
```

# Vraag 10 Test: Uitvoering
Op basis van de simulaties en de variantie-test, zien we dat de boxplot die wij bekwamen niet zo uitzonderlijk is, onder de hypothese van gelijke varianties. Daarom zullen we er voor kiezen om een t-test met gepoolde steekproefvariantie te gebruiken, omdat die in dat geval een iets hogere power heeft dan de t-test met Welch-benadering. Sla de test op in `t_test`

```{r}
t_test <-...
t_test
```

#Interpretatie

We vinden onvoldoende bewijs om te kunnen stellen dat het gemiddelde gewichtsverlies in de verbannen groep verschilt van het gemiddelde gewichtsverlies in de groep die niet verbannen werden op het 5% significantieniveau (p-waarde = $0.16$). We vinden een gemiddeld verschil in gewichtsverlies tussen de verbannen en de niet verbannen groep van 25.4g met een 95% BI van $-10.64$ tot $61.44$.

We kunnen op basis van deze observationele studie opnieuw geen causaal verband leggen tussen gewichtsverlies en het al dan niet verbannen worden uit de groep. Het is bijvoorbeeld mogelijk dat de individuen die verbannen worden reeds fysiek zwakker zijn dan de individuen die kunnen blijven, waardoor ze sowieso al een hogere kans zouden hebben om aan minder voedsel te raken en dus gewichtsverlies zouden lijden. Men zou een causaal verband kunnen leggen indien men willekeurig stokstaartjes in de groep zou laten blijven en stokstaartjes zouden forceren om de groep te verlaten.