Assignment3.Rmd

---
title: "Assignment 3 - Causal inference"
author: "RF"
date: "2/5/2020"
output: html_document
---

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

## Assignment 3 - Exploring causal inference issues

In this assignment we explore some issues related to multiple regressions (regressions with more than one predictor), and inferred (causal) relations between variables. N.B. the data is simulated (to make sure I know the actual mechanism generating it), but it's based on a real study. So bear with a longish introduction to get into the details of what we are doing and why it is important.

### Altercentric intrusion in schizophrenia

People with schizophrenia often report altered control and distinction of self-other representations: intrusive thoughts, hearing of voices, delusions of mind reading, paranoia, etc (a substantial portion of the psychotic symptoms experienced in schizophrenia). These have been variously attributed to hypermentalizing (over attribution of mental states to others), social impairment (over preoccupation with own thought processes), hyper socialization (inability to inhibit information from others), etc.

The current study investigates 1) whether schizophrenia is indeed related to altered control and distinction of self-other representations, in particular altercentric intrusions (inability to inhibit social information), and 2) whether these are related to the relevant psychotic symptoms. N.B. the actual study also investigates egocentric intrusion, do check the papers below if interested.

The task is a slightly modified version of this: https://www.ncbi.nlm.nih.gov/pubmed/20731512 You look at a picture with some dots visible to you, as well as with a different person with a different set of dots visible to them. The number of dots you see and that the other sees can be the same (congruent condition) or not (incongruent condition). You are tasked to indicate whether a given number (e.g. 3) matches the number of dots you see (and the dots visible to the other person are irrelevant to the task).


The tasks investigates altercentric intrusion: will your reaction time change according to whether the other person is seeing the same amount of dots as you, or not? The idea is that if you correctly inhibit social information, your reaction time should not change, as the information about the other person is not relevant. On the contrary, if you nevertheless use task irrelevant social information, you'll be slower at indicating whether 3 is the right number of dots when the other person sees a different amount of dots than you (conflicting information).
The bigger the difference between RTs in the congruent and incongruent condition the bigger the altercentric intrusion effect.

For each participant you have 6 variables: 1) ID, 2) AltercentricIntrusion (continuous score), 3) Diagnosis (schizophrenia vs. control), 4) VoiceHearing (severity of voice hearing symptoms, continuous score of the severity of the symptom as measured by a clinician), 5) MindReading (severity of delusions of mind reading, continuous score of the severity of the symptom as measured by a clinician); 6) Apathy (severity of lack of motivation in taking care of oneself, from washing to showing up at work, continuous score of the severity of the symptom as measured by a clinician).

The research questions you have to answer are the following:

## First part

Q1.1) Does schizophrenia involve altercentric intrusion? Define model and priors. Test the implications of your priors (prior predictive checks) and if needed adjust them. Run the model. Test the quality of the fitted model (posterior predictive checks). Assess the evidence in favor of an increased altercentric intrusion in schizophrenia. Report the model and the results, including plots.

# Q1

Based on the summary of the model, it seems that there is a difference in altercentric intrusion of 0.37 points between controls and patients, with patients scoring highest (as hypothesized). The model estimates the mean score of controls to be 3.86, and the sigma (the expected average error of the mean) to be 0.92.
( b = , CIs = , )
On average, controls show an effect of xxx (CIs = ) in Altercentric intrusion. Schizophrenic patients show an effect of xx (CIs = )
[Add plots of the effects]

Plotting the data against the model reveals that the model does not capture the large variance in both of the groups.

```{r}
data3 <- read.csv("Ass3.csv")
data3$Diagnosis <- plyr::revalue(as.character(data3$Diagnosis), 
                             c("0"="Controls", "1"="Schizophrenia"))
data3 <- data3 %>%
  mutate(
    ID = as.factor(ID),
    Diagnosis = as.factor(Diagnosis)
  )
library(pacman)
pacman::p_load(tidyverse, brms, rethinking, patchwork)
```


```{r}
# Define model
m1_f <- bf(AltercentricIntrusion ~ 0 + Diagnosis)


# Ask which priors need to be defined
get_prior(m1_f, data = data3, family = gaussian)

# Define priors
prior <- c(
  prior(normal(4, 1), class = b),
  prior(normal(1, 2), class = sigma)
)

# Draw the consequences of priors AND likelihood (why AND?)
m1_prior <- brm(m1_f,data3, family = gaussian, prior = prior, sample_prior= "only") # "only" means we are not running it on the data, but on samples from the prior

# Test implications of priors (prior predictive checks)
pp_check(m1_prior, nsamples = 100)

# Run the model
m1 <- brm(m1_f,data3,prior = prior, sample_prior= T)
summary(m1)
plot(m1)

# Quality test of the fitted model (posterior predictive checks)
pp_check(m1)

# plots, plots, plots
data3$Diagnosis <- as.factor(data3$Diagnosis)
  

data3 %>%
  ggplot(aes(x = Diagnosis, y = AltercentricIntrusion)) +
  geom_abline(intercept = fixef(m1)[1], 
              slope     = fixef(m1)[2]) +
  geom_point(shape = 1, size = 2, color = "pink") +
  theme_bw() +
  theme(panel.grid = element_blank()) +
  geom_jitter()


# Hypothesis testing
plot(hypothesis(m1,
                "DiagnosisSchizophrenia > DiagnosisControls"))
hypothesis(m1,
           "DiagnosisSchizophrenia > DiagnosisControls")
conditional_effects(m1)
plot(conditional_effects(m1), points=T)

```



Q1.2) Is altercentric intrusion related to specific symptoms *in the patients*? Identify which of the symptoms could be relevant. Should you include more than one symptom? Build models, priors, predictive checks. Assess the evidence and report models and results, including plots. Discuss whether the results make sense.

```{r}

# filter patients
data_p <- filter(data3, Diagnosis == "Schizophrenia")

# scale the data
data_p <- data_p %>% 
  mutate(
    AltercentricIntrusion = scale(AltercentricIntrusion),
    VoiceHearing = scale(VoiceHearing),
    MindReading = scale(MindReading),
    Apathy = scale(Apathy)
    )

# Define models
m2_f <- bf(AltercentricIntrusion ~ 1 + VoiceHearing)
m3_f <- bf(AltercentricIntrusion ~ 1 + MindReading)
m4_f <- bf(AltercentricIntrusion ~ 1 + Apathy)

# Ask which priors need to be defined
get_prior(m2_f, data = data_p, family = gaussian)

# Define priors
prior2 <-c(
  prior(normal(0, 1), class = Intercept), # Altercentric Intrusion. Intercept should be 0 because we scaled
  prior(normal(0, .3), class = b),  # Voice hearing
  prior(normal(1, 1), class = sigma))
prior3 <-c(
  prior(normal(1, 2), class = Intercept), # Altercentric Intrusion
  prior(normal(2, 1), class = b),  # Mind reading
  prior(normal(1, 1), class = sigma))
prior4 <-c(
  prior(normal(4, 1), class = Intercept), # Altercentric Intrusion
  prior(normal(2, 0.5), class = b),  # Apathy
  prior(normal(1, 1), class = sigma))


# Draw the consequences of priors AND likelihood (why AND?)
m2_prior <- brm(m2_f, data_p, prior = prior2, family = gaussian, sample_prior = "only")
m3_prior <- brm(m3_f, data_p, prior = prior3)
m4_prior <- brm(m4_f, data_p, prior = prior4)

# Test implications of priors (prior predictive checks)
pp_check(m2_prior, nsamples = 100)
pp_check(m3_prior, nsamples = 100)
pp_check(m4_prior, nsamples = 100)

## Run the models ##
# voice hearing model
m2 <- brm(
  formula = m2_f,
  data = data_p, 
  family = gaussian,
  prior = prior2, 
  sample_prior= T) 
summary(m2)
plot(m2)

 # mind reading model
m3 <- brm(
  formula= m3_f, 
  data= data_p, 
  family = gaussian,
  prior = prior2, 
  sample_prior = T)

# apathy model
m4 <- brm(
  formula= m4_f, 
  data = data_p, 
  family = gaussian,
  prior = prior2, 
  sample_prior = T)

plot(data_p$VoiceHearing ~ data_p$MindReading)

# Quality test of the fitted model (posterior predictive checks)
pp_check(m2, nsamples = 100)
pp_check(m3, nsamples = 100)
pp_check(m4, nsamples = 100)

# Check models for warnings
m2
m3
m4

## Hypothesis testing
# Voice hearing
plot(hypothesis(m2,
                "VoiceHearing > 0"))
hypothesis(m2,
           "VoiceHearing > 0")
conditional_effects(m2)
plot(conditional_effects(m2), points=T)

# Mind reading
plot(hypothesis(m3,
                "MindReading > 0"))
hypothesis(m3,
           "MindReading > 0")
conditional_effects(m3)
plot(conditional_effects(m3), points=T)

# Apathy
plot(hypothesis(m4,
                "Apathy > 0")) # negative --> the more apathetic, the less altercentric intrusion
hypothesis(m4,
           "Apathy < 0")
conditional_effects(m4)
plot(conditional_effects(m4), points=T) # the result may be credible (it is), but it doesn't explain a lot of the variance in the data


## Multiple regression ##

# Mind reading + Voice hearing
m5 <- 
  brm(data = data_p, 
      family = gaussian,
      AltercentricIntrusion ~ 1 + VoiceHearing + MindReading,
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)

print(m5)
mcmc_plot(m5)

# mind reading + voice hearing + apathy
m6 <- 
  brm(data = data_p, 
      family = gaussian,
      AltercentricIntrusion ~ 1 + VoiceHearing + MindReading + Apathy,
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)

print(m6)
mcmc_plot(m6) 


# hypothesis testing
plot(hypothesis(m5, "VoiceHearing > 0"))
hypothesis(m5, "VoiceHearing > 0")

plot(hypothesis(m5, "MindReading > 0"))
hypothesis(m5, "MindReading > 0")

# Once we have a multiple regression model, we change the question we're asking in the hypothesis testing. Before: if we know VH, does that tell us something about AI? Here: if we already know A and MR, does also knowing VH add something? (that is what "VoiceHearing > 0" means in this case). So the answer we get is, that if we already know A and MR, VH doesnt credibly add any information

## MODEL COMPARISON ## (maybe we're just overfitting)
# asking: which model best predicts AI without overfitting the data
m2 <- add_criterion(m2, criterion = "loo")
m3 <- add_criterion(m3, criterion = "loo")
m4 <- add_criterion(m4, criterion = "loo")
m5 <- add_criterion(m5, criterion = "loo")

# which model is estimate to have the lowest PE?
loo_compare(m2, 
            m3, 
            m4, 
            m5) # m4 (Apathy) seems to be the best

# model weights: if we assume that these are the only possible models and one of them is true, what is the probability of any of the models to be the true model?
loo_model_weights(m2, 
            m3, 
            m4, 
            m5) # apathy is the true model, if these are the only possible models (maybe we don't have the true model)


```

Given our domain knowledge ... we expect VH and MR (but not apathy) to be related to AI. 
Models predictin AI from single symptoms du not support these hypotheses
report: VH, MR, A

A model comparison approach indicates that the model predicting AI from A is the best model, minimizing estimated out-of-sample-error (stacking weight of 1). Adding other symptoms to the model with A does not improve generalizability of the model (stacking weights of 0)

The results do not support our hypotheses and would require a rethinking of the theoretical assumptions. 

## Second part

Q2.1) However, we know that the diagnosis is based on symptom assessment: if the overall sum of symptoms is severe enough, the participant gets a diagnosis. In other words, by selecting the patients, and including the symptoms in the model we might have inadvertently introduced an issue in our inference. Do try to draw a causal graph (Directed Acyclical Graph) of the variables and compare it with the types of causal graphs presented in the slides. Discuss which biases you might have introduced.


Q2.2.) Redesign your analysis following the graph and report how the results change

```{r}

# scale the data
d <- data3 %>% 
  mutate(
    AltercentricIntrusion = scale(AltercentricIntrusion),
    VoiceHearing = scale(VoiceHearing),
    MindReading = scale(MindReading),
    Apathy = scale(Apathy)
    )

# # # Define models # # #
# VOICE HEARING MODEL
m2.2 <- brm(
  data= d,
  family = gaussian,
  AltercentricIntrusion ~ 1 + VoiceHearing,
  prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)

# MIND READING MODEL
m3.2 <- brm(
  data= d,
  family = gaussian,
  AltercentricIntrusion ~ 1 + MindReading,
  prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)

# APATHY MODEL
m4.2 <- brm(
  data= d,
  family = gaussian,
  AltercentricIntrusion ~ 1 + Apathy,
  prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)



# Quality test of the fitted model (posterior predictive checks)
pp_check(m2.2, nsamples = 100)
pp_check(m3.2, nsamples = 100)
pp_check(m4.2, nsamples = 100)


## Hypothesis testing
# Voice hearing
plot(hypothesis(m2.2,
                "VoiceHearing > 0"))
hypothesis(m2.2,
           "VoiceHearing > 0")
m2.2
conditional_effects(m2.2)
plot(conditional_effects(m2.2), points=T)

# Mind reading
plot(hypothesis(m3.2,
                "MindReading > 0"))
hypothesis(m3.2,
           "MindReading > 0")
m3.2
conditional_effects(m3.2)
plot(conditional_effects(m3.2), points=T)

# Apathy
plot(hypothesis(m4.2,
                "Apathy > 0"))
hypothesis(m4.2,
           "Apathy > 0")
m4.2
conditional_effects(m4.2)
plot(conditional_effects(m4.2), points=T) 


## Multiple regression ##

# Mind reading + Voice hearing
m5.2 <- 
  brm(data = d, 
      family = gaussian,
      AltercentricIntrusion ~ 1 + VoiceHearing + MindReading,
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)

print(m5.2)
mcmc_plot(m5.2)

# mind reading + voice hearing + apathy
m6.2 <- 
  brm(data = d, 
      family = gaussian,
      AltercentricIntrusion ~ 1 + VoiceHearing + MindReading + Apathy,
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, .3), class = b),
                prior(normal(1, 2), class = sigma)),
      sample_prior = T)

print(m6.2)
mcmc_plot(m6.2) 


# hypothesis testing
plot(hypothesis(m5.2, "VoiceHearing > 0"))
hypothesis(m5.2, "VoiceHearing > 0")

plot(hypothesis(m5.2, "MindReading > 0"))
hypothesis(m5.2, "MindReading > 0")


## MODEL COMPARISON ## (maybe we're just overfitting)
# asking: which model best predicts AI without overfitting the data
m2.2 <- add_criterion(m2.2, criterion = "loo")
m3.2 <- add_criterion(m3.2, criterion = "loo")
m4.2 <- add_criterion(m4.2, criterion = "loo")
m5.2 <- add_criterion(m5.2, criterion = "loo")
m6.2 <- add_criterion(m6.2, criterion = "loo")

# which model is estimate to have the lowest PE?
loo_compare(m2.2, 
            m3.2, 
            m4.2, 
            m5.2,
            m6.2) # m5.2 (VH + MR) has the lowest PE

# model weights
loo_model_weights(m2.2, 
            m3.2, 
            m4.2, 
            m5.2,
            m6.2) # m5.2 is most probable to be the true model (stacking weights = .75)

```


## Third part

These issues are very difficult to think through, and not knowing the causal mechanisms generating the data in advance makes our inferences even more unreliable. To explore these issues, I recommend using simulations. In other words, defining a "true" model, generating data from it and assessing what different analyses would lead you to infer (and therefore which biases they might introduce). You can find the code I used to simulate your data below.

Q3.1) Look through the code and identify whether the results you have match the underlying truth. Discuss what you have learned.

Q3.2) OPTIONAL: is this a general pattern? Try varying the parameters (e.g. correlation values) and assess whether the new dataset(s) leads to the same biases in your analysis.



```{r}
pacman::p_load(MASS, tidyverse, psych)

seed <- 1981 # Defining a seed so the results are always the same
n <- 300 # Defining the amount of participants

SymptomCorr <- .2 # Defining the correlation of symptoms (as they tend to co-occur)
EffectCorrRel <- .2 # Defining the correlation between relevant symptoms and effect (Some symptoms are positively correlated with the effect)
EffectCorrIrrel <- 0 # Defining the correlation between irrelevant symptoms and effect (none)

# Creating the variance-covariance matrix for the variables we want to generate (3 symptoms, 1 effect)
Sigma <- matrix(data=c(1,SymptomCorr,SymptomCorr,EffectCorrRel,
                       SymptomCorr,1,SymptomCorr,EffectCorrRel,
                       SymptomCorr,SymptomCorr,1,EffectCorrIrrel,
                       EffectCorrRel,EffectCorrRel,EffectCorrIrrel,1),
                       nrow=4,ncol=4)

## Generate data from a multivariate (mvr) normal (n) distribution
d <- mvrnorm(n = n, # number of participant
        mu = c(1.2, 1.2, 1.2, 4), # mean of each variable
        Sigma) # variance co-variance matrix

# Giving meaningful names to variables and add ID
d <- data.frame(
  VoiceHearing = d[,1], 
  MindReading =  d[,2],
  Apathy =  d[,3], 
  AltercentricIntrusion = d[,4],
  ID = seq(nrow(d)))

# Assessing whether the participant has schizophrenia (high enough sum of symptoms)
# Here we choose participants scoring above 75% percentile (the most severe ones)
d$Diagnosis <- 0
d$Diagnosis[(d$VoiceHearing + d$MindReading + d$Apathy) > 
              quantile(d$VoiceHearing + d$MindReading + d$Apathy, .75)] <-1

## Plotting the relation between variables in schizophrenia
d1 <- d %>% subset(Diagnosis==1) %>% dplyr::select(-Diagnosis, -ID)
pairs.panels(d1)

## Plotting the relation between variables all participants
pairs.panels(dplyr::select(d,-Diagnosis, -ID))

write_csv(d, "data/Ass3.csv")
```