---
title: "Upotreba jackknife procedure za analizu osetljivosti rezultata na outliere"
output: 
  html_document:
    toc: true
    toc_float: 
      collapsed: false
    toc_depth: 3
    code_folding: hide
    theme: simplex
    highlight: tango
    df_print: kable
date: "`r Sys.Date()`"
---

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

Za ilustraciju jackknife procedure iskoristićemo eksperiment o uticaju socijalnog isključivanja na perifernu temperaturu. Originalna studija IJzerman et al. (2012) (https://doi.org/10.1016/j.actpsy.2012.05.002) je pokazala da telesna temperatura ispitanika (merena na prstu) opada u eksperimentalnoj situaciji socijalnog isključivanja, dok kod grupe koja je uključena nema promene u temperaturi. Trenutno je u planu velika replikaciona studija koja će pokušati da potvrdi ili opovrgne ovaj efekat u većem broju nezavisnih laboratorija (https://osf.io/v5s6d/). Za potrebe nove studije kreirali smo skript koji sadrži konfirmatorne analize, uključujući i jackknife postupak. Ovde je izdvojen jedan deo tog skripta.


# Pripremne aktivnosti

## Učitavanje potrebnih paketa

paket faux služi generisanju baza podataka koje imaju višeslojnu strukturu   
paket dyplyr služi lakšoj manipulaciji podacima  
paket lme4 služi modelovanju mešovitih modela  
paket pbkrtest služi poređenju fita mešovitih modela 

```{r loading packages, message=FALSE, warning=FALSE}

library(faux)
library(dplyr)
library(lme4)
library(pbkrtest)

```

## Generisanje simuliranog uzorka koji odgovara potrebnoj analizi i računanje svih potrebnih skorova

Kreiraćemo mali uzorak da bi se analize izvele relativno brzo.  
Na početku postavljamo seed da bismo uvek dobili iste rezultate kada pokrenemo skript.

```{r jackknife sample}

set.seed(2345)

d <- add_random(lab = 3) %>%
  add_random(ID = round(runif(3, min=10, max=20)), .nested_in = "lab") %>%
  add_random(point = 15) %>%
  add_between("ID", condition = c("inclusion","exclusion")) %>%
  add_recode("condition", "Condition", inclusion = 0, exclusion = 1) %>%
  add_ranef("ID", temp_baseline = 0.5)  %>%
  mutate(temp_baseline = 35 + temp_baseline) %>%
  add_ranef("point", temp_actual = 0.5, t_slope = 0.5, .cors = 0) %>%
  add_ranef(sigma = 1) %>%
  mutate(temp_actual = 35 + temp_actual + (0 + t_slope) * Condition + sigma) %>%
  add_between("point", tik = 0:14) %>%
  add_ranef(c("ID","point"), t = 1) %>%
  mutate(t = as.numeric(tik) * 10 + t)



d$temp <- d$temp_actual - d$temp_baseline # deviation of actual temperature from a person's baseline temperature

d$time <- round(d$t/10 - 0.49, 0) # rounding time to multiple of 10 (fully elapsed) seconds
d <- mutate(d, time = ifelse(time < 0, 0, time)) # transforming negative values to 0

d$condition_ref_zero <- d$Condition - 0.5 # centering the condition variable, i.e. zero as the reference group
d$condition_ref_included <- d$Condition # using the included group as the reference group
d$condition_ref_excluded <- d$Condition - 1 # using the excluded group as the reference group

d$time_ref_0 <- d$time - 18 # centering the time variable around the midpoint (180s)
d$time_ref_max <- d$time - 36 # centering the time variable around it's maximum value (360s)


```


# Glavni model

U pitanju je generalni mešani model (general mixed model) gde je periferna temperatura na kažiprstu zavisna varijabla, a prediktori su vreme merenja (kontinuirani prediktor), eksperimentalna manipulacija (isključenost vs. uključnenost) i njihova interakcija.
Kako je struktura podataka višeslojna (multilevel) - mere temperature su ugnježdene u ispitanike, pored fiksnih faktora, u model uključujemo i slučajne faktore. Jedan slučajni faktor odnosi se na intercept za ispitanike (nemaju svi istu prosečnu temperaturu), a drugi na efekat vremena na temperaturu (kod različitih ispitanika promene u temperaturi mogu biti različite brzine). Na kraju, kako postoji još jedan nivo ugnježdavanja - ispitanici su ugnježdeni u laboratorije - uključujemo i slučajan efekat laboratorije u model.


Glavna istraživačka hipoteza tiče se efekta interakcije vremena i eksperimentalne grupe. Kako bismo procenili njenu značajnost, pravimo dva modela - glavni istraživački model i model za poređenje. Model za poređenje je isti kao i istraživački, osim što ne sadrži interkaciju. Poređenjem fita ova dva modela dobijamo procenu veličine i značajnosti efekta.


Istraživački model

```{r confirmatory model}

m1<- lmer(temp ~  (time | ID) + time + condition_ref_zero + time:condition_ref_zero + (1 | lab), data=d)
summary(m1)

```

Model za poređenje i poređenje dva modela

```{r confirmatory-comparisson model}

m0<- lmer(temp ~  (time | ID) + time + condition_ref_zero + (1 | lab), data=d)
KRmodcomp(m1, m0)

```

Pošto smo generisali jako mali uzorak (za kompleksnu strukturu podataka) i slučajan - efekat nije značajan, ali to zanemarujemo.


# Analiza osetljivosti na outliere - jackknife

Da bismo proverili da li su rezultati osetljivi na outliere koristimo jackknife.  


Prvo, pravimo funkciju koja (iz prethodnog outputa) izvlači i prikazuje samo značajnost interakcije.

## p za interakciju
```{r jackknife model p}

get_p <- function(d) {
m1<- lmer(temp ~  (time | ID) + time + condition_ref_zero + time:condition_ref_zero + (1 | lab), data=d)
m0<- lmer(temp ~  (time | ID) + time + condition_ref_zero + (1 | lab), data=d)
k<-KRmodcomp(m1, m0)
p=k$stats['p.value']
p
}

```


Zatim, kreiramo funkciju koja isključuje samo po jednog ispitanika iz uzorka u datom trenutku, za sve ispitanike. Drugim rečima, kreiramo N uzoraka čija je veličina N-1 (i svi su različiti jer se svaki put isključuje drugi ispitanik iz uzorka) - jackknife postupak.

## jackknife postupak
```{r jackknife remove one}

unique(d$ID) # original sample - unique participant IDs

# jackknifed subsamples - unique participant IDs
for (i in unique(d$ID)) {  
dd <- d[which(d$ID!=i),]
print(unique(dd$ID))
}

```


Na kraju, koristimo jackknife postupak da procenimo da li su rezultati osetljivi na outliere tako što poredimo originalni rezultat sa rezultatima dobijenim na jackknife poduzorcima.


## primena jackknife postupka
```{r jackknife interaction, warning=FALSE, message=FALSE}

# p-value for the original sample
res <- cbind(out='none', pvalue=get_p(d))
print(res)

# p-values for the jackknifed subsamples
for (i in unique(d$ID)) {  
  dd<-d[which(d$ID!=i),]
  string <- paste('case ',i)
  res <- cbind(out=string, pvalue=get_p(dd))
  print(res)
}

```