R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

Uncontrols: Shock response in right BST

library(brms) # for the analysis

Loading required package: Rcpp
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Loading 'brms' package (version 2.12.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').

Attaching package: 㤼㸱brms㤼㸲

The following object is masked from 㤼㸱package:stats㤼㸲:

    ar

library(haven) # to load the SPSS .sav file
library(tidyverse) # needed for data manipulation.

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.3
v tidyr   1.0.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(RColorBrewer) # needed for some extra colours in one of the graphs
library(ggmcmc)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

library(ggthemes)
library(ggridges)

Model

\[Y \sim N(\mu, \sigma_{\epsilon}^{2})\] \[\mu = A + \alpha*T + \beta*S+ \epsilon\] Where TM, TD, SM, and SD are covariates.

T: Trait score
S: State score

Priors

\[ A \sim N(0,100)\] \[\alpha \sim N(0,100)\] \[\beta \sim N(0,100)\] \[\sigma_{\epsilon} \sim Cauchy(0,100)\]

setwd("C:/Users/Chirag/Box/Box/UMD/Project_UMD/eCON/RBA/uncon_rBNST_with_covariates")

df <- read.table('uncon_rBNST_with_covariates.txt',header=TRUE)
prior1 <- c(prior(normal(0,100),class=Intercept),
            prior(normal(0,100),class=b, coef="TRAIT"),
            prior(normal(0,100),class=b, coef="STATE"),
            prior(cauchy(0,100),class=sigma)
           )

bmod1 <- brm(Y ~ TRAIT + STATE, 
             data = df, 
             family = gaussian(),
             prior = prior1,
             warmup = 2000, iter = 5000,
             chains = 4,
             cores  = 2)

Compiling the C++ model
Start sampling

Model Summary

summary(bmod1)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Y ~ TRAIT + STATE 
   Data: df (Number of observations: 61) 
Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
         total post-warmup samples = 12000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     0.42      0.05     0.33     0.52 1.00    11378     8426
TRAIT        -0.06      0.06    -0.17     0.05 1.00    10321     8629
STATE         0.08      0.06    -0.03     0.19 1.00    10617     8382

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.38      0.04     0.32     0.46 1.00    10536     8696

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Posteriors

data <- fixef(bmod1,summary = FALSE)
write.csv(data,file = 'posteriors.csv')

# IMPORT LIBRARIES:
library(data.table)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.12.8 using 2 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attaching package: 㤼㸱data.table㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    between, first, last

The following object is masked from 㤼㸱package:purrr㤼㸲:

    transpose

library(ggplot2)
library(ggridges)
library(dplyr)
library(tidyr)
library(scales)

Attaching package: 㤼㸱scales㤼㸲

The following object is masked from 㤼㸱package:purrr㤼㸲:

    discard

The following object is masked from 㤼㸱package:readr㤼㸲:

    col_factor

library(viridis)

Loading required package: viridisLite

Attaching package: 㤼㸱viridis㤼㸲

The following object is masked from 㤼㸱package:scales㤼㸲:

    viridis_pal

rois <- colnames(data)
#data$X <- NULL

nobj=dim(data)[1]

# rename columns with ROI list


data_stats = data.frame(1:length(rois))

# create ROI column instead of numerics to match threat table above

data_stats$ROI <- rois
data_stats$mean <- colMeans(data)
data_stats$P <- colSums(data > 0)/nobj
data_stats$Pn <- data_stats$P


for (i in 1:length(rois)) {
  if (data_stats$P[i]<.5){data_stats$Pn[i]=1-data_stats$P[i]}
}


# this will order the distributions correctly
data_stats <- data_stats[order(data_stats$mean),]



data_trans <- as.data.frame(t(as.matrix(data)))
data_trans <- tibble::rownames_to_column(data_trans, "ROI")
data_trans$X <- 1:nrow(data_trans)




# merge values & stats into one table by ROI

data_merge <- merge(data_stats, data_trans, by = "ROI")
data_merge <- data_merge[order(data_merge$X),]

# Transform data into long form


# Melt dataframe by ROI
library(data.table)

data_long <- melt(data_trans, id=c("ROI","X"))

The melt generic in data.table has been passed a data.frame and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(data_trans). In the next version, this warning will become an error.

#threat_long <- melt(threat_trans, id=c("ROI","X.x"))


data_long <- data_long[order(data_long$X),]


# This is incredibly clunky, but for the sake of time adding stats by ensuring orders are all the same and repeating each value 2000 times... Tried a few different methods with no success for some reason. 

data_long$mean <- rep(data_merge$mean, each = nobj)
data_long$P <- rep(data_merge$P, each =nobj)
data_long$Pn <- rep(data_merge$Pn, each =nobj)

# set your labels here so you don't have to change within the plot below: 
format(round(1, 2), nsmall = 2)

[1] "1.00"

y.axis.labs <- format(round(data_stats$P,2), nsmall = 2)                        # y axis labels
sec.y.axis.labs <- data_stats$ROI                             # second y axis labels (probabilities)


 # X AXIS LABELS NEED TO CHANGE TO CORRESPOND TO DATA SET! UNCOMMENT WHICHEVER MATCHES

# Uncomment for THREAT
x.axis.labs <- c( "Shock < 0", "0", "Shock > 0")                 # x axis labels  THREAT
x.labs.pos <- c(-0.2, 0, 0.2)                                             # x axis position THREAT


# Uncomment for VALENCE
#x.axis.labs <- c("Neutral > Positive", "0", "Positive > Neutral")         # x axis labels  VALENCE
#x.labs.pos  <- c(-0.1,0,0.2)                                              # x axis position VALENCE

# Uncomment for INTERACTION
#x.axis.labs <- NULL                                                       # x axis labels  INTERACTION, not sure what to put.
#x.labs.pos <- NULL                                                        # a axis position INTERACTION, change when labels decided

#graph.title <- data.name                                   # graph title 
legend.title <- "Probability"                              # legend title
y.axis.title <- NULL                                       # for now ...
x.axis.title <- NULL                                       # for now...

# GRAPH DATA 
dataset <- data_long
x.values <- data_long$value                               # x values
y.values <- data_long$ROI                                 # y values
y.values.RO <- data_long$value                            # values to reorder Y by
distrib.fill <- data_long$Pn                              # fill graph with probabilities
group <- data_long$ROI

# Other aspects
gradient.colors <- c("yellow","#C9182B","#41245C")        # change gradient colors here (currently, yellow - purple)
label.size <- 12                                         # adjust ROI and probability y-axis font size
title.size <- 28                                          # adjust graph title size 
x.axis.size <- 12                                        # adjust x-axis label sizes

# here is where you can change information about the graph and add other characteristics using ggplot and ggridges


ggplot(dataset, aes(x = x.values, y = as.numeric(reorder(y.values, y.values.RO)), 
                    fill = distrib.fill, group = group)) +                        # scale = spacing, alpha = transparency
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = 2, 
                      alpha = .95, 
                      scale = 2.55,
                      color = "black",
                      size = .35
                     ) +            # divide into two quantiles (show mean)
  geom_vline(xintercept = 0, linetype="solid",color = "black",alpha = .95, size = .45) +    #create line at X = 0
  scale_fill_gradientn(colors = viridis_pal(direction = 1, option = "inferno")(20),                         # set gradient
                       limits = c(.85,1),                                # which probabilites matter?
                       na.value = "#909497",                                # if not in limits, gray out
                       name = legend.title) +                            # name legend
  scale_y_continuous(breaks = 1:length(rois),                            # A VERY HACK-Y WAY TO HAVE TWO Y AXES W DISCRETE DATA
                     expand = c(0,0),
                     labels = y.axis.labs,                               # Trick ggplot into thinking data is continuous...
                     sec.axis = sec_axis(~.,                             # Second axis to show probabilities
                                         breaks = 1:length(rois),
                                         labels = sec.y.axis.labs)) +
  #theme_ridges(font_size = label.size, grid = TRUE, center_axis_labels = TRUE) +  # theme info
  #ggtitle(graph.title)+ 
  #theme_bw() +# graph title
 #theme_ridges(grid = FALSE) + 
  theme(   
    panel.background = element_blank(),
    legend.position = "none",
    #panel.grid.major.y = element_line(color = "grey"), 
    plot.title = element_text(hjust = 0.5, size = title.size),            # plot title size and position
    axis.text.y = element_text(size=label.size),                               # y-axis text size
    axis.line.x = element_line(color = "gray"),
    axis.text.y.right = element_text(size = label.size),                  # y-axis info for right axis
    axis.text.x = element_text(size = x.axis.size),
    #plot.margin = unit(c(0,0,0,0), "cm"),
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    legend.title.align = 5)+
  guides(shape = guide_legend(label.position = "bottom",
                              title.positon = "bottom",
                              title.vjust = 0.4)) +          
  labs(
    x = NULL,                                                         # Add or not add X and Y labels
    y = NULL) +
  scale_x_continuous(breaks = x.labs.pos, labels = c(x.axis.labs))

#ggsave(file = paste0(strsplit(type,'.RData'),".jpg", sep=""), width=4.5, height=3, units = "in", dpi = 900)

# preview will look smooshed, be sure to look at the saved .jpg before making any changes to text size, etc!
#graph