DecisionAnalysis.Rmd

---
title: "Decision analysis"
output:
  html_document:
    df_print: paged
---

Analyze decisions across teams.

```{r setup, echo=FALSE,message=FALSE}

use_checkpoint = TRUE

s = Sys.info()
if (s['sysname'] == 'Darwin'){
  use_checkpoint = FALSE
}

if (use_checkpoint) {
  library(checkpoint)
  checkpointDir <- '/checkpoint'
  checkpoint("2019-08-13", checkpointLocation = checkpointDir)
  }


library(plyr)
library(tidyverse)
library(lmerTest)
library(lme4)
library(emmeans) 
library(pscl)
library(arm)
library(MuMIn)
library(multcomp)
library(multcompView)
library(optimx)
library(xtable)
library(psych)
library(knitr)
library(sessioninfo)
library(boot)


basedir = Sys.getenv('NARPS_BASEDIR')
if (basedir == ""){
   # use default
   basedir = "/data"
}

# get details about installed packages
package_info <- devtools::package_info()
write.table(package_info, 
            file=paste(basedir,                     
                       'metadata/R_package_info.txt',
                       sep='/'))
```

## Data setup 

Load the data and clean up some variables.  Requires metadata that is created by PrepareMetadata.ipynb

```{r loadData}
# load and clean up data
cat(sprintf('using basedir: %s', basedir))

narps_df <- read_csv(paste(basedir,
                           'metadata/all_metadata.csv',
                           sep='/'))

narps_df <- narps_df %>% 
  mutate(Confidence = as.ordered(narps_df$Confidence))

narps_df$testing[is.na(narps_df$testing)] <- "Other"
narps_df$teamID = as.factor(narps_df$teamID)
narps_df$varnum = as.factor(narps_df$varnum)
```

Estimate correlation between reported smoothing kernel and estimated image smoothness.

```{r smoothCorr}
cor(narps_df$fwhm,narps_df$smoothing_coef,
    use='pairwise.complete',
    method='spearman')
```

## Make decision summary table

Compute data for Table 1 in paper

```{r}
desc <- c('Positive parametric effect of gains in the vmPFC (equal indifference group)',
'Positive parametric effect of gains in the vmPFC (equal range group)',
'Positive parametric effect of gains in the ventral striatum (equal indifference group)',
'Positive parametric effect of gains in the ventral striatum (equal range group)',
'Negative parametric effect of losses in the vmPFC (equal indifference group)',
'Negative parametric effect of losses in the vmPFC (equal range group)',
'Positive parametric effect of losses in the amygdala (equal indifference group)',
'Positive parametric effect of losses in the amygdala (equal range group)',
'Greater positive response to losses in amygdala for equal range group vs. equal indifference group')

decision_df <- narps_df %>%
  dplyr::select(Decision, varnum, Similar, Confidence) 

decision_summary <- decision_df %>%
  group_by(varnum) %>%
  summarise(pDecision = mean(Decision), 
            medianSimilar = median(Similar),
            madSimilar = mad(Similar, constant=1),
            medianConfidence = median(as.numeric(Confidence)),
            madConfidence = mad(as.numeric(Confidence), constant=1))
decision_summary$Description <- desc
decision_summary <- decision_summary[c('varnum',
                                       'Description',
                                       'pDecision',
                                       'medianSimilar',
                                       'madSimilar',
                                       'medianConfidence',
                                       'madConfidence')]
write.table(decision_summary,
            file=paste(basedir,"figures/Table1.tsv",sep='/'),
            sep='\t')

decision_summary

```

## Check independent variables

Note that the correlations are not interpretable for the factorial variables - the plot is shown simply to provide some insight into the relationships between the different variables and their distributions.

```{r vis_model}

df_for_vis <- narps_df %>% 
  drop_na(fwhm,used_fmriprep_data,package,testing) %>%
  dplyr::select(c(fwhm,used_fmriprep_data,package,testing,smoothing_coef, movement_modeling))

pairs.panels(df_for_vis)

```

## Descriptive stats


```{r dataSetup}
hyp_df = narps_df %>% drop_na(fwhm,
                              used_fmriprep_data,
                              package,
                              testing)
hyp_df <- hyp_df %>%  mutate(
    package = factor(package, c('SPM', 'FSL', 'AFNI', 'Other')),
    testing = factor(testing),
    used_fmriprep_data = factor(used_fmriprep_data)) %>%
  mutate(testing = recode_factor(hyp_df$testing, ARI = 'Other', randomise = 'nonparametric', permutations = 'nonparametric'))

save(hyp_df,file=paste(basedir,'output/hyp_df.RData', sep='/'))
```

### FWHM

```{r}
print(mean(narps_df$fwhm,na.rm=TRUE))
print(min(narps_df$fwhm,na.rm=TRUE))
print(max(narps_df$fwhm,na.rm=TRUE))

```

### fMRIPrep

```{r tableFmriprep}

table(hyp_df$used_fmriprep_data)/9

```

### Package


```{r tablePackage}

table(hyp_df$package)/9

```

### Testing

```{r tableTesting}

table(hyp_df$testing)/9

```
### Movement

```{r tableMovement}

table(hyp_df$movement_modeling)/9

```

## Models

First run mixed model across full dataset to assess overall effects on hypothesis acceptance. We use nlminb as the optimizer because the model failed to converge using the standard optimizer in lmer.

```{r fullModel}

m_hyp_full = glmer(Decision ~ varnum + fwhm + used_fmriprep_data + package + testing + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_full_model <- summary(m_hyp_full)
print(summary_full_model)
sink(paste(basedir,'figures/decision_model_summary.txt', sep='/'))
print(summary_full_model)
sink() 
save(m_hyp_full, summary_full_model,file=paste(basedir,'output/full_model_summary.RData', sep='/'))

r.squaredGLMM(m_hyp_full)

```

### Nonparametric bootstrap
The bootMer package performs resampling for lmer models, but it uses a parametric bootstrap which we did not feel adequately addressed concerns about non-normality.  Thus, we implemented a nonparametric bootstrap, in which we resample over teams.

First we generate wide data, for resampling.

```{r}
decisions_wide <- hyp_df %>%
  mutate(Hypothesis = str_c("Hyp", varnum, sep="")) %>%
  dplyr::select(teamID, Decision, Hypothesis) %>%
  spread(Hypothesis,Decision)

team_info <- hyp_df %>% 
  filter(varnum==1) %>%
  dplyr::select(fwhm, used_fmriprep_data,package,
                testing, movement_modeling,teamID)

decisions_wide <- join(decisions_wide,
                       team_info,
                       by='teamID')
```

Function to run analysis within bootstrap.  For the regressors that are parametric or binary, we can simply compute a bootstrap confidence interval on their parameter estimate.  For the two variables that are factorial, we can't do this because some bootstrap samples will not contain all of the same values for the factors, and thus the models are not equivalent.  Instead, for those variables we perform a model comparison and save the model comparison statistics so that we can compute a bootstrap estimate on those (ala https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5482523/).  

```{r message=FALSE,warning=FALSE}
bs_fun <- function(df, indices){
  df <- df[indices,]
  control = glmerControl(
              optimizer ='optimx', 
              optCtrl=list(method='nlminb'))
  # we take in a wide data frame, but must
  # convert it back to long form for lmer
  df_long <- df %>% gather(key='Hypothesis',
              value='Decision',
              -movement_modeling,-teamID,-fwhm,
              -used_fmriprep_data,-package,-testing)
  
  model_full <- glmer(Decision ~ fwhm + 
                        used_fmriprep_data + 
                        movement_modeling + 
                        Hypothesis + 
                        package + 
                        testing  + (1|teamID),
            data = df_long,family=binomial,
            control = control)
  
  # get coefs for all but factorial effects
  fixed_effect_coefs = fixef(model_full)[2:4]
  
   # do model comparison for Hypothesis effect
  model_nohyp <- glmer(Decision ~ fwhm + 
                             used_fmriprep_data + 
                             movement_modeling + 
                             package + 
                             testing  + (1|teamID),
            data = df_long,family=binomial,
            control = control)

  # perform model comparison and save statistics
  anova_nohyp <- anova(model_nohyp,model_full)
  fixed_effect_coefs['hyp_BICdiff'] = 
    diff(anova_nohyp$BIC)
  fixed_effect_coefs['hyp_AICdiff'] = 
    diff(anova_nohyp$AIC)

  # do model comparison for package effect
  model_nopackage <- glmer(Decision ~ Hypothesis + 
                             fwhm + 
                             used_fmriprep_data + 
                             movement_modeling + 
                             testing  + (1|teamID),
            data = df_long,family=binomial,
            control = control)

  # perform model comparison and save statistics
  anova_nopackage <- anova(model_nopackage,model_full)
  fixed_effect_coefs['package_BICdiff'] = 
    diff(anova_nopackage$BIC)
   fixed_effect_coefs['package_AICdiff'] = 
    diff(anova_nopackage$AIC)
 
  # do model comparison for testing effect
   model_notesting <- glmer(Decision ~ Hypothesis + 
                              fwhm + 
                              used_fmriprep_data + 
                              movement_modeling + 
                              package +  (1|teamID),
            data = df_long,family=binomial,
            control = control)
    
  # perform model comparison and save statistics
  anova_notesting <- anova(model_notesting,model_full)
  fixed_effect_coefs['testing_BICdiff'] = 
    diff(anova_notesting$BIC)
  fixed_effect_coefs['testing_AICdiff'] = 
    diff(anova_notesting$AIC)

fixed_effect_coefs
}

```


```{r}
# perform bootstrap sampling with 1000 replications
Out <- boot(data=decisions_wide, statistic=bs_fun, R=1000)

# save bootstrap samples
save(Out,file=paste(basedir,                     
                     'output/npboot_output.RData',
                     sep='/'))

```

Get confidence intervals, first for betas for parametric regressors.

```{r}
# set up data frame to save results
rownames = c('fwhm',
            'used_fmriprep_data',
            'movement_modeling', 
            'Hypothesis (deltaBIC)',
            'Hypothesis (deltaAIC)',
            'package (deltaBIC)',
            'package (deltaAIC)',
            'testing (deltaBIC)',
            'testing (deltaAIC)') 


percentile_np_ci <- data.frame(lower=array(NA, 
                                    length(rownames)),
                               mean=array(NA, 
                                    length(rownames)),
                               upper=array(NA,
                                    length(rownames)),
                               exceedence=array(NA,
                                    length(rownames)),
                               row.names=rownames)

# get means for parameters across amples
npboot_means = apply(Out$t,2,mean)

# loop over parametric regressors and get percentile CI
for (i in 1:3){
  ci_raw <- boot::boot.ci(boot.out = Out, type = c("perc"), index = i)
  percentile_np_ci[i,] <- c(ci_raw$percent[4],npboot_means[i],ci_raw$percent[5],NA)
}

```

Get CI for model comparisons for factorial regressors (package and testing).

```{r}
# get ci for BIC diff

for (i in 4:9){
  ci_raw <- boot::boot.ci(boot.out = Out, type = c("perc"), index = i)
  # use lower p because we want to know how often full model
  # is better than reduced model
  percentile_np_ci[i,] <- c(ci_raw$percent[4],npboot_means[i],ci_raw$percent[5],mean(Out$t[,i]<0))
}

kable(percentile_np_ci)
write.table(percentile_np_ci,
            file=paste(basedir,"output/nonparametric_bootstrap_CI.tsv",sep='/'),
            sep='\t')

```


Get odds ratios for factorial variables.

```{r}
# from https://stackoverflow.com/questions/26417005/odds-ratio-and-confidence-intervals-from-glmer-output
odds_ratios <- exp(fixef(m_hyp_full))
# this uses the immensely faster but less accurate Wald method for confidence intervals
cc <- confint(m_hyp_full,parm="beta_",method="Wald")  
ctab <- cbind(est=odds_ratios,exp(cc))
kable(ctab)
write.table(ctab,
            file=paste(basedir,"figures/OddsRatios.tsv",sep='/'),
            sep='\t')

```


### Model comparisons

Estimate a set of models leaving out each variable of interest, so that we can then use model comparison to estimate the effect sizes.  This is particularly necessary for the factor variables since they are not represented by a single variable in the model.

#### Model without hypothesis

```{r hypModel}
m_hyp_nohyp = glmer(Decision ~ fwhm + used_fmriprep_data + package + testing + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))
summary(m_hyp_nohyp)
anova_hyp<- anova(m_hyp_full,
                  m_hyp_nohyp,
                  test='Chisq')
print(anova_hyp)

# compute delta r-squared between this model and full model
delta_r2_hyp <- r.squaredGLMM(m_hyp_full) - r.squaredGLMM(m_hyp_nohyp) 
print(delta_r2_hyp)
```

#### Model without smoothing

```{r nosmoothModel}
m_hyp_nosmooth = glmer(Decision ~ varnum + used_fmriprep_data + package + testing + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))
summary(m_hyp_nosmooth)
anova_smoothing <- anova(m_hyp_full,
                         m_hyp_nosmooth,
                         test='Chisq')
print(anova_smoothing)

emtrends(m_hyp_full,specs='fwhm',var='fwhm')

# compute delta r-squared between this model and full model
delta_r2_smoothing <- r.squaredGLMM(m_hyp_full) - r.squaredGLMM(m_hyp_nosmooth) 
print(delta_r2_smoothing)
```

#### Model without fmriprep

```{r noprepModel}
m_hyp_noprep = glmer(Decision ~ varnum + fwhm + package + testing + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))
summary(m_hyp_noprep)
anova_fmriprep <- anova(m_hyp_full,
                        m_hyp_noprep,
                        test='Chisq')
print(anova_fmriprep)

emmeans(m_hyp_full,'used_fmriprep_data')

# compute delta r-squared between this model and full model
delta_r2_fmriprep <- r.squaredGLMM(m_hyp_full) - r.squaredGLMM(m_hyp_noprep) 

print(delta_r2_fmriprep)
```

#### Model without software package

```{r nopackageModel}


m_hyp_nopackage = glmer(Decision ~ varnum + used_fmriprep_data + fwhm + testing + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb'))) 
summary(m_hyp_nopackage)
anova_package <- anova(m_hyp_full,
                       m_hyp_nopackage,
                       test='Chisq')
print(anova_package)

leastsquare = emmeans(m_hyp_full,
                      'package')
multcomp::cld(leastsquare, 
    level=.05)

# compute delta r-squared between this model and full model
delta_r2_package <- r.squaredGLMM(m_hyp_full) - r.squaredGLMM(m_hyp_nopackage)
print(delta_r2_package)
```

#### Model without testing method

```{r notestingModel}
  
m_hyp_notesting = glmer(Decision ~ varnum + used_fmriprep_data + fwhm + package + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))   
summary(m_hyp_notesting)
anova_testing <- anova(m_hyp_full,
                       m_hyp_notesting,test='Chisq')
print(anova_testing)

leastsquare = emmeans(m_hyp_full,
                      'testing')
multcomp::cld(leastsquare, 
    level=.05)

# compute delta r-squared between this model and full model
delta_r2_testing <- r.squaredGLMM(m_hyp_full) - r.squaredGLMM(m_hyp_notesting)
print(delta_r2_testing)
```

#### Model without movement modeling

```{r nomovementModel}
  
m_hyp_nomovement = glmer(Decision ~ varnum + used_fmriprep_data + fwhm + package  + testing + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))   
summary(m_hyp_nomovement)
anova_movement <- anova(m_hyp_full,
                       m_hyp_nomovement,test='Chisq')
print(anova_movement)

# compute delta r-squared between this model and full model
delta_r2_movement <- r.squaredGLMM(m_hyp_full) - r.squaredGLMM(m_hyp_nomovement)
print(delta_r2_movement)
```

#### Create table for paper with modeling results

```{r}
summary_df = data.frame(Effects=c('Hypothesis','Estimated smoothness','Used fMRIPprep data','Software package','Multiple correction method','Movement modeling'))
summary_df['Chi-squared'] = c(anova_hyp$Chisq[2],
                              anova_smoothing$Chisq[2],
                              anova_fmriprep$Chisq[2],
                              anova_package$Chisq[2],
                              anova_testing$Chisq[2],
                              anova_movement$Chisq[2])
summary_df['P value (parametric)'] = c(anova_hyp$"Pr(>Chisq)"[2],
                              anova_smoothing$"Pr(>Chisq)"[2],
                              anova_fmriprep$"Pr(>Chisq)"[2],
                             anova_package$"Pr(>Chisq)"[2],
                              anova_testing$"Pr(>Chisq)"[2],
                              anova_movement$"Pr(>Chisq)"[2])
summary_df['Bootstrap CI'] = c(NA,
                              sprintf('(%0.2f - %0.2f)',
                                      percentile_np_ci['fwhm','lower'],
                                      percentile_np_ci['fwhm','upper']),
                              sprintf('(%0.2f - %0.2f)',
                                      percentile_np_ci['used_fmriprep_data','lower'],
                                      percentile_np_ci['used_fmriprep_data','upper']),
                              NA,
                              NA,
                              sprintf('(%0.2f - %0.2f)',
                                      percentile_np_ci['movement_modeling','lower'],
                                      percentile_np_ci['movement_modeling','upper']))
summary_df['Model selection probability (BIC)'] = c(percentile_np_ci['Hypothesis (deltaBIC)','exceedence'],
                              NA,
                              NA,
                              percentile_np_ci['package (deltaBIC)','exceedence'],
                              percentile_np_ci['testing (deltaBIC)','exceedence'],
                              NA)
summary_df['Model selection probability (AIC)'] = c(percentile_np_ci['Hypothesis (deltaAIC)','exceedence'],
                              NA,
                              NA,
                              percentile_np_ci['package (deltaAIC)','exceedence'],
                              percentile_np_ci['testing (deltaAIC)','exceedence'],
                              NA)
summary_df['Delta R^2'] = c(delta_r2_hyp[2,1],
                            delta_r2_smoothing[2,1],
                            delta_r2_fmriprep[2,1],
                            delta_r2_package[2,1],
                            delta_r2_testing[2,1],
                            delta_r2_movement[2,1])
summary_df <- summary_df %>%
  mutate(`Chi-squared` = formatC(summary_df$'Chi-squared',
                                 digits=2,format='f'),
    `Delta R^2` = formatC(summary_df$'Delta R^2',
                        digits=2,format='f')
    )
summary_df
write.table(summary_df,
            file=paste(basedir,
                       "figures/ModelingSummaryTable.tsv",sep='/'),
            quote=FALSE,sep='\t', row.names = FALSE)
```

#### Modeling excluding Hypotheses 7-9

Due to the low amount of variability in outcomes for Hypotheses 7-9, we assessed whether excluding those hypotheses would affect modeling outcomes.

```{r excludeModel}

hyp_df_exc <- hyp_df %>%
  filter(varnum %in% c(1,2,3,4,5,6))

m_hyp_full_exc = glmer(Decision ~ varnum + fwhm + used_fmriprep_data + package + testing + movement_modeling + (1|teamID),
            data = hyp_df_exc,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_full_model <- summary(m_hyp_full_exc)
print(summary_full_model)


```


#### Applied vs. estimated smoothness

Full model using specified smoothing filter rather than estimated smoothness.

### FWHM

```{r}
print(median(narps_df$smoothing_coef,na.rm=TRUE))
print(min(narps_df$smoothing_coef,na.rm=TRUE))
print(max(narps_df$smoothing_coef,na.rm=TRUE))

```

```{r kernelModel}
m_hyp_full_tsc = glmer(Decision ~ varnum + smoothing_coef + used_fmriprep_data + package + testing + movement_modeling + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))
summary(m_hyp_full_tsc)

```


### Separate analyses for each variable of interest

#### Hypothesis 

```{r varnum_only}
m_hyp_varnum = glmer(Decision ~ varnum + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_varnum <- summary(m_hyp_varnum)
print(summary_varnum)

```

#### Estimated smoothness

```{r fwhm_only}
m_hyp_fwhm = glmer(Decision ~ fwhm + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_fwhm <- summary(m_hyp_fwhm)
print(summary_fwhm)

```

#### Applied smoothing kernel

```{r kernel_only}
m_hyp_kernel = glmer(Decision ~ smoothing_coef + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_kernel <- summary(m_hyp_kernel)
print(summary_kernel)

```

#### fMRIprep

```{r fmriprep_only}
m_hyp_fmriprep = glmer(Decision ~ used_fmriprep_data + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_fmriprep <- summary(m_hyp_fmriprep)
print(summary_fmriprep)

```

#### Package

```{r package_only}
m_hyp_package = glmer(Decision ~ package  + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_package <- summary(m_hyp_package)
print(summary_package)

```

#### Testing

```{r testing_only}
m_hyp_testing = glmer(Decision ~ testing + (1|teamID),
            data = hyp_df,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

summary_testing <- summary(m_hyp_testing)
print(summary_testing)

```

### Separate analyses for each hypothesis

For each hypothesis, we ask whether decisions are associated with fwhm and use of fmriprep.  None of these survive Bonferroni correction.

```{r runModelsFunction}

runModels = function(hyp, data_df){
  pvals = c()

  m_hyp_full = bayesglm(Decision ~ fwhm + used_fmriprep_data + package + testing + movement_modeling,
                   data = data_df,family=binomial)
  
  m_hyp_nosmooth = bayesglm(Decision ~ used_fmriprep_data + package + testing,
                   data = data_df,family=binomial)
  #cat('testing effect of smoothing\n')
  a = anova(m_hyp_full,
            m_hyp_nosmooth,
            test='Chisq')
  pvals = c(pvals,unlist(a)['Pr(>Chi)2'])
  
  m_hyp_noprep = bayesglm(Decision ~ fwhm + package + testing,
                   data = data_df,family=binomial)
  #cat('testing effect of fmriprep\n')
  a=anova(m_hyp_full,m_hyp_noprep,test='Chisq')
  pvals = c(pvals,unlist(a)['Pr(>Chi)2'])
  
  m_hyp_nopackage = bayesglm(Decision ~ used_fmriprep_data + fwhm + testing,
                   data = data_df,family=binomial)
  #cat('testing effect of package\n')
  a=anova(m_hyp_full,
          m_hyp_nopackage,
          test='Chisq')
  pvals = c(pvals,unlist(a)['Pr(>Chi)2'])
  
  
  m_hyp_notesting = bayesglm(Decision ~ used_fmriprep_data + fwhm + package,
                   data = data_df,family=binomial)
  #cat('testing effect of testing\n')
  a=anova(m_hyp_full,
          m_hyp_notesting,
          test='Chisq')
  pvals = c(pvals,unlist(a)['Pr(>Chi)2'])
  return(pvals)
}
```


```{r runAllModels}
all_pvals = c()
for (hyp in 1:9){
  pv = runModels(hyp, hyp_df)
  all_pvals = rbind(all_pvals,pv)
}

```

### Assess consistency across packages

```{r}
var_by_package <- hyp_df %>% group_by(varnum, package) %>% summarize(varDecision = var(Decision))
var_model <- lm(varDecision ~ varnum + package, data = var_by_package)
var_model_nopackage <- lm(varDecision ~ varnum, data = var_by_package)
summary(var_model)
anova(var_model, var_model_nopackage)
emmeans(var_model,'package')

```

### Assess relation between smoothness and movement modeling

```{r smoothnessMovement}

smoothness_full = lmer(fwhm ~ varnum +  movement_modeling +  smoothing_coef + (1|teamID),
            data = hyp_df)

summary_smoothness <- summary(smoothness_full)
print(summary_smoothness)

```

### Examine effects of various factors on pattern distance

```{r loadPatternDistance}

patterndist_df = read_csv(paste(basedir,
          'metadata/median_pattern_corr.csv',
          sep='/'),skip=0)
names(patterndist_df)=c('teamID','median_corr')

# use mean FWHM across hypotheses

mean_fwhm <- narps_df %>% 
  group_by(teamID) %>% 
  summarize(fwhm = mean(fwhm))


merged_df = join(patterndist_df,
                 narps_df %>% filter(varnum==1),
                 by=c('teamID')) %>%
  dplyr::select(-fwhm)
  
merged_df = join(merged_df, mean_fwhm, by=c('teamID'))
```

### fit model

```{r r2z, echo=FALSE}
r2z = function(r){
    # fisher transform
    z=0.5*log((1.0+r)/(1.0-r))
    z[is.na(z)]=0
    return(z)
}
```

```{r fitPatternDistanceModel}

merged_df <- merged_df %>%
  mutate(z_median_corr = r2z(median_corr))

merged_df = merged_df %>% 
  drop_na(fwhm,used_fmriprep_data,package,testing)

dist_hyp_full = lm(z_median_corr ~ fwhm + used_fmriprep_data + package + testing + movement_modeling,
            data = merged_df)

s <- summary(dist_hyp_full)
print(s)

```

```{r fitnoSmoothDistanceModel}
dist_hyp_nosmooth <- lm(z_median_corr ~ used_fmriprep_data + package + testing + movement_modeling,
            data = merged_df)

anova(dist_hyp_full,
      dist_hyp_nosmooth,
      test='Chisq')

emtrends(dist_hyp_full,specs='fwhm',var='fwhm')
s$r.squared - summary(dist_hyp_nosmooth)$r.squared

```

```{r fitNoPrepDistanceModel}
dist_hyp_noprep = lm(z_median_corr ~ fwhm + package + testing + movement_modeling,
            data = merged_df)

anova(dist_hyp_full,
      dist_hyp_noprep,
      test='Chisq')

emmeans(dist_hyp_full,'used_fmriprep_data')
s$r.squared - summary(dist_hyp_noprep)$r.squared

```

```{r fitNoPackageDistanceModel}
dist_hyp_nopackage = lm(z_median_corr ~ used_fmriprep_data + fwhm + testing + movement_modeling,
            data = merged_df) 

anova(dist_hyp_full,
      dist_hyp_nopackage,
      test='Chisq')

s$r.squared - summary(dist_hyp_nopackage)$r.squared
leastsquare = emmeans(dist_hyp_full,'package')
multcomp::cld(leastsquare, 
    level=.05)


```

```{r fitNoTestingDistanceModel}
dist_hyp_notesting = lm(z_median_corr ~ used_fmriprep_data + fwhm + package + movement_modeling,
            data = merged_df)   

anova(dist_hyp_full,
      dist_hyp_notesting,
      test='Chisq')

s$r.squared - summary(dist_hyp_notesting)$r.squared

leastsquare = emmeans(dist_hyp_full,
                      'testing')

multcomp::cld(leastsquare, 
    level=.05)


```

```{r fitNoMovementDistanceModel}
dist_hyp_nomovement = lm(z_median_corr ~ used_fmriprep_data + fwhm + package + testing,
            data = merged_df)   

anova(dist_hyp_full,
      dist_hyp_nomovement,
      test='Chisq')

s$r.squared - summary(dist_hyp_nomovement)$r.squared


```

## Similarity as function of fmriprep use


## Make supplementary table 5

```{r mkSuppTable5}

df_SuppTable4 <- data.frame(
  Effects = c('Hypothesis', 
              'Smoothness', 
              'Used fMRIprep data', 
              'Software package', 
              'Multiple testing correction')) %>%
  mutate(ChiSquared = NA,
         Pvalue = NA,
         delta_r2 = NA)

df_SuppTable4[1,2:4] = c(anova_hyp$`Chisq`[2],
                         anova_hyp$`Pr(>Chisq)`[2],
                         delta_r2_hyp[2,1]
                         )
df_SuppTable4[2,2:4] = c(anova_smoothing$`Chisq`[2],
                         anova_smoothing$`Pr(>Chisq)`[2],
                         delta_r2_smoothing[2,1]
                         )
df_SuppTable4[3,2:4] = c(anova_fmriprep$`Chisq`[2],
                         anova_fmriprep$`Pr(>Chisq)`[2],
                         delta_r2_fmriprep[2,1]
                         )
df_SuppTable4[4,2:4] = c(anova_package$`Chisq`[2],
                         anova_package$`Pr(>Chisq)`[2],
                         delta_r2_package[2,1]
                         )
df_SuppTable4[5,2:4] = c(anova_testing$`Chisq`[2],
                         anova_testing$`Pr(>Chisq)`[2],
                         delta_r2_testing[2,1]
                         )

kable(df_SuppTable4)

write.table(df_SuppTable4,
            file=paste(basedir,"figures/SuppTable4.tsv",sep='/'),
            sep='\t')



```


# create tidy version of pairwise map similarity

```{r makeSimdata}
simdata <- data.frame(varnum=integer(),
                      team1=character(),
                      team2=character(),
                      similarity = double())

for (hyp in c(1,2,5,6,7,8,9)){
  # load pairwise similarity
  sim_df <- read_csv(
    paste(basedir,
    'output/correlation_unthresh',
    sprintf('spearman_unthresh_hyp%d.csv',hyp),
    sep='/'))
  rownames(sim_df) <- sim_df$X1
  sim_df$X1 <- NULL
  for (i in 1:dim(sim_df)[1]){
    for (j in 1:dim(sim_df)[1]){
      if (j > i){
        team1_fmriprep = hyp_df %>% dplyr::filter(teamID==rownames(sim_df)[i]) %>% slice(1) %>% pull(used_fmriprep_data)
        team2_fmriprep = hyp_df %>% dplyr::filter(teamID==names(sim_df)[j]) %>% slice(1) %>% pull(used_fmriprep_data)
        
        newdata <- data.frame(
          varnum=hyp,
          team1=rownames(sim_df)[i],
          team2=names(sim_df)[j],
          similarity = as.numeric(sim_df[i, j]),
          fmriprep = sum(c(team1_fmriprep=='Yes', team2_fmriprep=='Yes')))
        simdata <- rbind(simdata,newdata)
      }
    }
  }
  
}
simdata <- simdata %>%
  mutate(similarity_z = r2z(simdata$similarity),
         fmriprep = as.factor(simdata$fmriprep))

z2r <- function(z) {
  return((exp(2.0*z) - 1)/(exp(2.0*z) + 1))
}

simdata %>% group_by(fmriprep) %>% summarise(m = z2r(mean(similarity)))
```


Performing inference on the differnces between these groups would require implementing a permutation test approach.  Given how small the numeric differences are between groups, I decided to forego this analysis for now.

### Decision vs confidence/similarity

How do confidence and similarity ratings relate to hypothesis decisions?

```{r}
# load data
narps_df_allteams <- read_csv(paste(basedir,
                           'metadata/narps_metadata_all_teams.csv',
                           sep='/'))
narps_df_allteams$Hypothesis <- as.factor(narps_df_allteams$varnum)

m_hyp_conf_sim = glmer(Decision ~ Hypothesis + Similar + Confidence + (1|teamID),
            data = narps_df_allteams,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))
summary(m_hyp_conf_sim)

m_hyp_conf_nosim = glmer(Decision ~ Hypothesis + Confidence + (1|teamID),
            data = narps_df_allteams,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

delta_r2_sim <- r.squaredGLMM(m_hyp_conf_sim) - r.squaredGLMM(m_hyp_conf_nosim)
print(sprintf('Delta r^2 for similarity: %0.3f', delta_r2_sim[2,1]))

m_hyp_noconf_sim = glmer(Decision ~ Hypothesis + Similar + (1|teamID),
            data = narps_df_allteams,family=binomial,
            control = glmerControl(optimizer ='optimx', 
                                   optCtrl=list(method='nlminb')))

delta_r2_conf <- r.squaredGLMM(m_hyp_conf_sim) - r.squaredGLMM(m_hyp_noconf_sim)
print(sprintf('Delta r^2 for confidence: %0.3f', delta_r2_conf[2,1]))

group_means <- narps_df_allteams %>% group_by(Decision) %>% summarize(mean_confidence=mean(Confidence),mean_similarity=mean(Similar))
kable(group_means)
```

Compute the Spearman correlation of the distance of the outcome (the proportion of teams that reported a significant result for each hypothesis) from 0.5 (i.e., how consistent the results were across teams) and the mean confidence rating across hypotheses.

```{r}
hyp_summary_df <- narps_df_allteams %>%
  group_by(Hypothesis) %>%
  summarize(mean_distance = abs(0.5 - mean(Decision)),
            mean_confidence = mean(Confidence),
            mean_similarity = mean(Similar))

cor.test(hyp_summary_df$mean_confidence,hyp_summary_df$mean_distance,method='spearman')

cor.test(hyp_summary_df$mean_similarity,hyp_summary_df$mean_distance,method='spearman')

```