Starting an EDA on the mult_wave part of the citizen_shield project and looking into a bunch of networks linking the determinants with the least amount of missing data

2022-06-28 20:58:43 +02:00
parent ac27be6bab
commit d3d4571467
13 changed files with 5723 additions and 1 deletions
@@ -0,0 +1,403 @@
+---
+title: "unified SEM and Multilevel VAR Estimation all determinants"
+output: html_notebook
+---
+
+# Outline
+This script covers the following steps
+
+A. N = 1 model - unified Structural Equation Model (uSEM)
+
+- Set up the uSEM model
+
+- Check model summary (model fit statistics)
+
+- Data visualization and interpretation of results
+
+B. N = all Multilevel VAR Model
+
+- Set up the mlVAR model
+
+- Check model summary(model fit statistics)
+
+- Data visualization and interpretation of results
+
+C. Selected Readings
+
+- for uSEM model: Yang, X., Ram, N., Gest, S., Lydon, D., Conroy, D. E., Pincus, A. L., & Molenaar, P. C. M. (2018). Socioemotional dynamics of emotion regulation and depressive symptoms: A person-specific network approach. Complexity, 2018, Article ID 5094179. doi: 10.1155/2018/5094179 [Open Access https://www.hindawi.com/journals/complexity/2018/5094179/]
+
+- for mlVAR model: Bringmann LF, Vissers N, Wichers M, Geschwind N, Kuppens P, Peeters F, et al. (2013) A Network Approach to Psychopathology: New Insights into Clinical Longitudinal Data. PLoS ONE 8(4): e60188. https://doi.org/10.1371/journal.pone.0060188 [Open Access https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0060188]
+
+### Preliminaries
+#### Loading Libraries
+
+Loading libraries used in this script
+
+```{r}
+# Check to see if necessary packages are installed, and install if not
+packages <- c("psych", "pompom", "mlVAR", "bootnet", "psychonetrics", "GGMncv", "EstimateGroupNetwork", "relaimpo")
+if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
+  install.packages(setdiff(packages, rownames(installed.packages())))  
+}
+# Load packages
+library(tidyr)
+library(lubridate)
+library(stringr)
+library(psych) #for data description
+library(plyr) #for data manipulation
+library(dplyr)
+library(ggplot2) #for data visualization
+# library(pompom) #for uSEM
+# library(mlVAR) #for mlVAR models
+library(bootnet)
+library(psychonetrics)
+library(EstimateGroupNetwork)
+```
+
+### Loading Data
+
+```{r}
+fs_demo_columns <- c('id', 'demographic_age', 'demographic_age_factor',
+       'demographic_education', 'demographic_gender',
+       'demographic_gender_factor', 'demographic_income',
+       'demographic_living_with', 'demographic_region',
+       # 'demographic_underage_children',
+       'demographic_underage_children_factor', 'fsd_end', 
+       # 'fsd_id',
+       # 'fsd_no', 
+       'fsd_round', 
+       'fsd_start', 'fsd_vnk', 
+       # 'fsd_vr',
+       'fsd_weight')
+
+```
+
+```{r}
+df <- read.csv("../citizen_shield/data/kp_df_eng_ordinal.csv", header=TRUE, row.names = "X")
+remove_na <- function(DF, n=0) {
+  DF[, colSums(is.na(DF)) <= n]
+}
+
+df <- remove_na(df, n=4000)
+
+```
+
+```{r}
+# describeBy(df, group="user_id")
+# tmp <- df %>%
+#   group_by(user_id) %>%
+#   summarise_at(vars(feature_list), funs(sd(., na.rm=TRUE)))
+```
+
+
+
+#### create a feature list based on the column names of the data frame
+
+```{r}
+feature_list <- colnames(df %>% select(-fs_demo_columns))
+```
+
+#### Show the head and description of the overall data frame
+
+```{r}
+head(df)
+```
+
+```{r}
+# describeBy(df, group = "user_id")
+# describe(df)
+```
+
+#### Show the head and descrption of an example user 
+
+```{r}
+example_round <- 1
+data_indiv <- df[df$fsd_round == example_round, ]
+# head(data_indiv)
+# describe(data_indiv)
+```
+
+
+```{r}
+plot_df <- data_indiv %>%
+  select(c(id, fsd_round), all_of(feature_list[1:20])) %>%
+  gather(key = "variable", value = "value", -c(id, fsd_round))
+```
+
+#### plotting intraindividual change 
+
+```{r}
+#plotting intraindividual change 
+ggplot(data = plot_df,
+       aes(x = id, y=value, group= fsd_round)) +
+  #first variable
+  geom_line(aes(color = variable)) + 
+  #plot layouts
+  scale_x_continuous(name="Arbitrary Time") +
+  scale_y_continuous(name="Raw Values") +  
+  theme_classic() +
+  theme(axis.title=element_text(size=14),
+        axis.text=element_text(size=14),
+        plot.title=element_text(size=14, hjust=.5)) +
+  ggtitle(example_round)
+```
+#### Normalize the data frame (consider min max maybe)
+
+```{r}
+# standardize specific data columns (not the id or time variables in first 3 columns)
+data_indiv[feature_list] <- lapply(data_indiv[feature_list], 
+                             function(x) c(scale(x, center=TRUE, scale=TRUE)))
+# describe(data_indiv)
+```
+
+#### plotting normalized intraindividual change 
+```{r}
+plot_df <- data_indiv %>%
+  select(c(id, fsd_round), all_of(feature_list[1:20])) %>%
+  gather(key = "variable", value = "value", -c(id, fsd_round))
+```
+
+
+```{r}
+#plotting intraindividual change 
+ggplot(data = plot_df,
+       aes(x = id, y=value, group= fsd_round)) +
+  #first variable
+  geom_line(aes(color = variable)) + 
+  #plot layouts
+  scale_x_continuous(name="Arbitrary Time") +
+  scale_y_continuous(name="Raw Values") +  
+  theme_classic() +
+  theme(axis.title=element_text(size=14),
+        axis.text=element_text(size=14),
+        plot.title=element_text(size=14, hjust=.5)) +
+  ggtitle(example_round)
+```
+
+
+Now we see that all the variables are in standardized form.
+
+### Check the data
+
+It is useful to check that there are data in all columns. If any one of the variables is all missing (or has no variance), the model cannot be fit. Missing data on a few observations within a column is ok.
+
+```{r}
+# check column missing
+na_col <- 0
+for (col in 1:ncol(data_indiv)) {
+  if (sum(is.na(data_indiv[,col])) == nrow(data_indiv)){
+        na_col <- na_col + 1
+  }
+}
+na_col
+```
+
+All columns are reported.
+
+### Describing the data.
+```{r}
+describe(df)
+```
+
+### Not standardizing the data
+
+For the following networks, the data are kept in their original form.
+
+```{r}
+## Impute the feature_list missingness
+# imp.cart <- mice::mice(df[, feature_list], method="cart", printFlag = FALSE)
+# df[, feature_list] <- mice::complete(imp.cart)
+complete_df <- df[complete.cases(df[, c("id", "fsd_round", feature_list)]), c("id", "fsd_round", feature_list)]
+```
+
+### ggmModSelect and EBICglasso networks
+There is no accounting for the hierarchical nature of the data here. Unregularized Gaussian Graphical Model ("ggmModSelect"; GGM) using the glasso algorithm and stepwise model selection. Gaussian Markov random field estimation using graphical LASSO and extended Bayesian information criterion ("EBICglasso") to select optimal regularization parameter.
+
+```{r}
+net_modSelect <- estimateNetwork(complete_df[feature_list], 
+              default = "ggmModSelect",
+              stepwise = FALSE,
+              corMethod = "cor")
+```
+
+```{r}
+net_thresh <- estimateNetwork(complete_df[feature_list],
+              tuning = 0, # EBICglasso sets tuning to 0.5 by default
+              default = "EBICglasso",
+              threshold = TRUE,
+              corMethod = "cor")
+```
+
+```{r}
+Layout <- qgraph::averageLayout(net_modSelect, net_thresh)
+layout(t(1:2))
+plot(net_modSelect, layout = Layout, title = "ggmModSelect", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+plot(net_thresh, layout = Layout, title = "Thresholded EBICglasso", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```
+
+Here the principal direction is forced - this rescales variables according to the sign of the first eigen-vector. This will lead to most correlations to be positive (positive manifold), leading to negative edges to be substantively interpretable. (not sure this is preferable in this instance as the variables are not all from the same questionnaire).
+
+```{r}
+net_modSelect_rescale <- estimateNetwork(complete_df[feature_list], 
+              default = "ggmModSelect",
+              stepwise = FALSE,
+              principalDirection = TRUE)
+net_thresh_rescale <- estimateNetwork(complete_df[feature_list],
+              tuning = 0, 
+              default = "EBICglasso",
+              threshold = TRUE,
+              principalDirection = TRUE)
+layout(t(1:2))
+plot(net_modSelect_rescale, layout = Layout, 
+     title = "ggmModSelect", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+plot(net_thresh_rescale, layout = Layout, 
+     title = "Thresholded EBICglasso", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```
+```{r}
+qgraph::centralityPlot(
+  list(
+    ggmModSelect = net_modSelect_rescale,
+    EBICGlasso_thresh = net_thresh_rescale
+  ), include = "ExpectedInfluence"
+)
+```
+
+```{r}
+boots <- bootnet(net_thresh_rescale, statistics = "ExpectedInfluence", 
+                 nBoots = 100, nCores = 2, type = "case")
+plot(boots, statistics = "ExpectedInfluence") + 
+  theme(legend.position = "none")
+```
+### relative importance network
+There is no accounting for the hierarchical nature of the data here
+
+```{r}
+# net_relimp <- estimateNetwork(complete_df[feature_list],
+#               default = "relimp",
+#               normalize = FALSE)
+# net_relimp2 <- estimateNetwork(complete_df[feature_list],
+#               default = "relimp",
+#               normalize = FALSE,
+#               structureDefault = "ggmModSelect",
+#               stepwise = FALSE # Sent to structureDefault function
+#     )
+```
+
+
+```{r}
+# Layout <- qgraph::averageLayout(net_relimp, net_relimp2)
+# layout(t(1:2))
+# plot(net_relimp, layout = Layout, title = "Saturated", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038")
+# plot(net_relimp2, layout = Layout, title = "Non-saturated", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038")
+```
+### graphicalVAR
+This approach is not really valid as it uses all time points as if they are from the same person which is not the case here
+
+```{r}
+# Estimate model:
+gvar <- estimateNetwork(
+  complete_df, default = "graphicalVAR", vars = feature_list,
+  tuning = 0, dayvar = "fsd_round", nLambda = 10
+)
+```
+
+
+```{r}
+Layout <- qgraph::averageLayout(gvar$graph$temporal,
+                                gvar$graph$contemporaneous)
+layout(t(1:2))
+plot(gvar, graph = "temporal", layout = Layout,
+     title = "Temporal", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+plot(gvar, graph = "contemporaneous", layout = Layout,
+     title = "Contemporaneous", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```
+
+```{r}
+# gvar_boot <- bootnet(gvar, nBoots = 2, nCores = 2)
+```
+
+```{r}
+# plot(gvar_boot, graph = "contemporaneous", plot = "interval")
+```
+
+### Mixed Graphical model
+
+```{r}
+net_mgm <- estimateNetwork(complete_df[feature_list],
+              default = "mgm",
+              type="g",
+              level=1
+              # type=c("g", "g", "g", "g", "c", "c", "c", "g"),
+              # level=c(1, 1, 1, 1, 12, 7, 15, 1)
+              )
+```
+
+```{r}
+plot(net_mgm, layout = Layout,
+     title = "mgm", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```
+
+```{r}
+par(mar=c(5,1,10,1)+7)
+for (checkpoint in unique(complete_df[, "fsd_round"])) {
+  net_mgm <- estimateNetwork(complete_df[complete_df["fsd_round"]==checkpoint, feature_list], 
+              default = "mgm",
+              type="g",
+              level=1)
+  plot(net_mgm, layout = Layout, 
+       # title = paste("mgm network", "checkpoint =", checkpoint), 
+       edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+  
+  title(main=paste("Mixed Graphical Model network", "checkpoint =", checkpoint), line=16.25)
+}
+```
+
+
+### Partial Correlation network
+
+```{r}
+net_par_cor <- estimateNetwork(complete_df[feature_list], 
+              default = "pcor")
+```
+
+```{r}
+plot(net_par_cor, layout = Layout, 
+     title = "Partial Correlation", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```
+
+```{r}
+par(mar=c(5,1,10,1)+7)
+for (checkpoint in unique(complete_df[, "fsd_round"])) {
+  net_par_cor <- estimateNetwork(complete_df[complete_df["fsd_round"]==checkpoint, feature_list], 
+              default = "pcor")
+  plot(net_par_cor, layout = Layout, 
+       edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+  
+  title(main=paste("Partial Correlation network", "checkpoint =", checkpoint), line=16.25)
+}
+```
+
+
+### Correlation network
+
+```{r}
+net_cor <- estimateNetwork(complete_df[feature_list], 
+              default = "cor")
+```
+
+```{r}
+plot(net_cor, layout = Layout, 
+     title = "Correlation", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```
+
+### GGMncv network
+
+```{r}
+net_GGMncv <- estimateNetwork(complete_df[feature_list], 
+              default = "GGMncv")
+```
+
+```{r}
+plot(net_GGMncv, layout = Layout, 
+     title = "GGMncv", edge.labels=TRUE, posCol="#306fbe", negCol="#e58038", label.scale.equal=TRUE, label.cex=10)
+```