########################################## ############### Models LMM ############### ########################################## rm(list=ls()) ### library ## library(psych) library(car) library(lme4) library(lmerTest) library(plyr) library(ggplot2) # load bdd setwd("/Users/cheval/switchdrive/Boris/Recherche/Recherche en cours/MT EEG/Zenodo/data management") load ("data_for_LMM.RData") ## transform variable as numeric y <- as.character(data_main_all_final$Stimuli) yr <- revalue(y, c("PA" = "-0.5", "SB" = "0.5")) data_main_all_final$stim_num <- as.numeric(yr) y <- as.character(data_main_all_final$Action) yr <- revalue(y, c("Approach" = "-0.5", "Avoid" = "0.5")) data_main_all_final$App_avoidance_num <- as.numeric(yr) #### Model with the IPAQ data_main_all_final$hab_MVPA_c <- data_main_all_final$MVPA - mean(data_main_all_final$MVPA, na.rm=T) data_main_all_final$hab_MVPA_Z <- data_main_all_final$hab_MVPA_c / sd(data_main_all_final$MVPA, na.rm=T) data_main_all_final$hab_MVPA_Z <- as.numeric(data_main_all_final$hab_MVPA_Z ) data_main_all_final$subject <- as.factor(data_main_all_final$subject) data_main_all_final$images <- as.factor(data_main_all_final$images) data_main_all_final$RT_rel_neutral R1_LMM_MVPA <- lmer(RT_rel_neutral ~ 1 + stim_num*App_avoidance_num*hab_MVPA_Z + (App_avoidance_num + stim_num | subject) + (1 | images), data=data_main_all_final, na.action = na.omit) summary(R1_LMM_MVPA) vcov(R1_LMM_MVPA) ################ ### diagnostics ############### # plot residual plot (R1_LMM_MVPA) # linearity in each variable ggplot(data.frame(x1=data_main_all_final$stim_num, pearson=residuals(R1_LMM_MVPA, type="pearson")), aes(x=x1, y =pearson)) + geom_point() + theme_bw() ggplot(data.frame(x1=data_main_all_final$App_avoidance_num, pearson=residuals(R1_LMM_MVPA, type="pearson")), aes(x=x1, y =pearson)) + geom_point() + theme_bw() ggplot(data.frame(x1=data_main_all_final$hab_MVPA_Z, pearson=residuals(R1_LMM_MVPA, type="pearson")), aes(x=x1, y =pearson)) + geom_point() + theme_bw() #### Multicollinearity # fonction vif pour mer (pris sur github, donc a verifier mais lauteur considere dit que # se calcul la multicolinearite pour lmer) # reste juste à adapter le nom de la fonction à la fin vif.mer <- function (fit) { ## adapted from rms::vif (variance inflation factor) v <- vcov(fit) nam <- names(fixef(fit)) ## exclude intercepts ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)")) if (ns > 0) { v <- v[-(1:ns), -(1:ns), drop = FALSE] nam <- nam[-(1:ns)] } d <- diag(v)^0.5 v <- diag(solve(v/(d %o% d))) names(v) <- nam v } vif.mer(R1_LMM_MVPA) # normality of residuals qqnorm(residuals(R1_LMM_MVPA)) hist(residuals(R1_LMM_MVPA)) # leverage ggplot(data.frame(lev=hatvalues(R1_LMM_MVPA), pearson =residuals(R1_LMM_MVPA, type="pearson")), aes(x=lev, y=pearson)) + geom_point () + theme_bw() #### simple slopes # -1SD data_main_all_final$hab_MVPA_Z_minus1SD <- data_main_all_final$hab_MVPA_Z + 1 R1_LMM_MVPA_minus1SD <- lmer(RT_rel_neutral ~ 1 + stim_num*App_avoidance_num*hab_MVPA_Z_minus1SD + (App_avoidance_num | subject) + (1 | images), data=data_main_all_final, na.action = na.omit) summary(R1_LMM_MVPA_minus1SD) # +1SD data_main_all_final$hab_MVPA_Z_plus1SD <- data_main_all_final$hab_MVPA_Z - 1 R1_LMM_MVPA_plus1SD <- lmer(RT_rel_neutral ~ 1 + stim_num*App_avoidance_num*hab_MVPA_Z_plus1SD + (App_avoidance_num | subject) + (1 | images), data=data_main_all_final, na.action = na.omit) summary(R1_LMM_MVPA_plus1SD)