library(arsenal)
library(car)
library(ggplot2)
library(GGally)
library(coin)
library(knitr)
library(markdown)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(visdat) #para evaluar datos faltantes
library(corrplot) #para graficar correlaciones
library(emmeans)
library(bestglm)
library(randomForest)
library(dplyr)
library(caret) #para k-fold cross validation

1 Objetivo principal del estudio

Analizar el funcionamiento ejecutivo en personas adultas cursando episodios depresivos melancólicos y no-melancólicos, y estudiar el efecto de otras variables clínicas y demográficas sobre el rendimiento.

2 Carga de la base de datos

La base de datos proviene de un estudio que recolecto información sobre personas adultas que cursan episodios depresivos melancólicos y no-melancólicos. Se recolectaron datos demográficos y clínicos de aspectos psiquiátricos y cognitivos de los pacientes. La base de datos completa contiene 58 variables.

#setwd("C:/Users/USER/Desktop/Taller/Trabajo final/TP final")
datos<-read.csv("base_taller_final")

Para el presente taller se seleccionaron alguna de ellas por tener relevancia para el objetivo del trabajo. Así, creamos un nuevo dataframe llamado melanco con las variables de interés más una variable sumaria para analizar el rendimiento cognitivo de los pacientes llamada FE (funciones ejecutivas) compuesta por el rendimiento promediado de tres pruebas ejecutivas: variables TMTB + FS + FF de la base original.

melanco <- datos[ , c(3:5, 7, 8, 11, 14, 15, 18, 20:28, 49:52,54:57)]
melanco$FE<-(melanco$TMTB +melanco$FS +melanco$FF)/3
melanco <- melanco[, c(1:10, 12:19, 27)]

3 Descripción de las variables

La base a analizar contiene 19 variables que se detallan a continuación:

VARIABLE RESPUESTA

FE: Continua. Desempeño estandarizado de los participantes en pruebas de evaluacion de funciones ejecutivas ajustado por edad y escolaridad. Mejor puntuación indica mejor rendimiento.

VARIABLE EXPLICATIVA DE MAYOR INTERES PARA EL ESTUDIO

Melanco_SMPI: Categorica. Diagnostico de melancolia segun criterios SMPI (1 = Melancolico; 0 = No-melancolico).

OTRAS VARIABLES EXPLICATIVAS Y/O DE CONTROL

Sexo: Categorica. 1 = Mujer; 0 = Hombre
Edad: Continua. Edad en años.
Educacion: Continua. Años de escolaridad formal.
TBP: Categorica. Tipo de trastorno (1 = TBP, 0 = TDM)
Edad_inicio: Continua. Edad (en años) en la que inicio la enfermedad.
Num_depresiones: Discreta. Cantidad de episodios depresivos.
Intentossuicidio: Discreta. Cantidad de intentos de suicidios.
Internaciones: Discreta. Cantidad de internaciones previas.
Psi_actuales: Categorica. Presencia de sintomas psicoticos en el EDM actual (1 = Si; 0 = No).
YMRS: Ordinal. Puntaje en la escala de mania.
MADRS: Ordinal. Puntaje en la escala de depresion.
IFD_B_BZD: Ordinal (5 niveles). Exposicion actual a benzodiacepinas medida en la escala IFD.
IFD_B_AD: Ordinal (5 niveles). Exposición actual a antidepresivos medida en la escala IFD.
IFD_B_EA: Ordinal (5 niveles). Exposición actual a estabilizadores del animo medida en la escala IFD.
IFD_B_AP: Ordinal (5 niveles). Exposición actual a antipsicóticos medida en la escala IFD.
IFD_B_Ach: Ordinal (2 niveles). Exposición actual a anticolinérgicos medida en la escala IFD.
TMTA: Continua. Rendimiento estandarizado en una prueba de evaluación atencional ajustado por edad y escolaridad. Puntajes más altos indican mejor rendimiento.

3.1 Puesta a punto de las variables para el analisis

Pasaje de las variables categóricas (dicotómicas y ordinales con pocos valores, por sugerencia de Adriana) a factores.

melanco$Sexo <- factor(melanco$Sexo, levels = c(0,1),
                       labels = c("Hombre", "Mujer"))
melanco$TBP <- factor(melanco$TBP, levels = c(0,1),
                      labels = c("TDM", "TBP"))
melanco$Melanco_SMPI <- factor(melanco$Melanco_SMPI, levels = c(0,1),
                               labels = c("No melancolico", "Melancolico"))

melanco$Psi_actuales<- factor(melanco$Psi_actuales, levels = c(0,1),
                                  labels = c("No", "Si"))
melanco$IFD_B_BZD<- factor(melanco$IFD_B_BZD)
melanco$IFD_B_EA<- factor(melanco$IFD_B_EA)
melanco$IFD_B_AD<- factor(melanco$IFD_B_AD)
melanco$IFD_B_AP<- factor(melanco$IFD_B_AP)
melanco$IFD_B_Ach<- factor(melanco$IFD_B_Ach)

4 Análisis exploratorio

4.1 Variables categoricas

La base contiene 9 variables categoricas: Sexo, TBP, Melanco_SMPI, Psi_actuales, IFD_B_BZD, IFD_B_AD, IFD_B_EA, IFD_B_AP, IFD_B_Ach. Se muestran gráficos exploratorios de la distribucion en función de la variable explicativa de interes Melanco_SMPI.

Figura 1.

4.2 Variables continuas o de conteo.

La base contiene 10 variables numericas (continuas o de conteo): Edad, Educacion, Edad_inicio, Num_depresiones, Intentossuicidios, Internaciones, YMRS, MADRS, TMTA, FE. Se muestran graficos exploratorios de la distribucion en funcion de la variable explicativa de interes Melanco_SMPI.

Figura 2.

4.3 Comparaciones bivariadas y tabla 1

Se compararan todas las variables con la explicativa principal del estudio Melanco_SMPI.

4.3.1 Variables categoricas

Evaluacion de supuestos para chi cuadrado

lapply(melanco[,c("Sexo", "TBP", "Psi_actuales", "IFD_B_BZD",
                  "IFD_B_AD", "IFD_B_EA", "IFD_B_AP", "IFD_B_Ach")],
       function(x) chisq.test(x, melanco$Melanco_SMPI)$expected)

## $Sexo
##         melanco$Melanco_SMPI
## x        No melancolico Melancolico
##   Hombre       26.66667    9.333333
##   Mujer        73.33333   25.666667
## 
## $TBP
##      melanco$Melanco_SMPI
## x     No melancolico Melancolico
##   TDM       51.85185    18.14815
##   TBP       48.14815    16.85185
## 
## $Psi_actuales
##     melanco$Melanco_SMPI
## x    No melancolico Melancolico
##   No       90.37037    31.62963
##   Si        9.62963     3.37037
## 
## $IFD_B_BZD
##    melanco$Melanco_SMPI
## x   No melancolico Melancolico
##   0      12.592593    4.407407
##   2       8.148148    2.851852
##   3       7.407407    2.592593
##   4      47.407407   16.592593
##   5      24.444444    8.555556
## 
## $IFD_B_AD
##    melanco$Melanco_SMPI
## x   No melancolico Melancolico
##   0      51.851852  18.1481481
##   2      14.074074   4.9259259
##   3      16.296296   5.7037037
##   4      15.555556   5.4444444
##   5       2.222222   0.7777778
## 
## $IFD_B_EA
##    melanco$Melanco_SMPI
## x   No melancolico Melancolico
##   0      54.814815   19.185185
##   2      28.888889   10.111111
##   3       8.148148    2.851852
##   4       5.185185    1.814815
##   5       2.962963    1.037037
## 
## $IFD_B_AP
##    melanco$Melanco_SMPI
## x   No melancolico Melancolico
##   0      10.370370    3.629630
##   2      60.000000   21.000000
##   3      16.296296    5.703704
##   4      10.370370    3.629630
##   5       2.962963    1.037037
## 
## $IFD_B_Ach
##    melanco$Melanco_SMPI
## x   No melancolico Melancolico
##   0      97.777778  34.2222222
##   2       2.222222   0.7777778

No se cumple el supuesto de 80% o más de las celdas con frecuencias absolutas esperadas mayores o iguales a 5 para sintomas psicóticos actuales (Psi_actuales) y todas las variables de la escala IFD (IFD_B_BZD, IFD_B_AD, IFD_B_EA, IFD_B_AP, IFD_B_Ach). Si se cumple para Sexo y TBP.

4.3.2 Variables continuas o de conteo

Evaluacion de supuestos para prueba t (normalidad e igualdad de varianzas).

Normalidad

Metodo analitico

lapply(melanco[,c("Edad", "Edad_inicio", "Educacion", "Num_depresiones",
                  "Intentossuicidio", "Internaciones", "YMRS", "MADRS",
                  "TMTA", "FE")],
       function(x) shapiro.test(lm(x ~ melanco$Melanco_SMPI)$residuals ))

## $Edad
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.97614, p-value = 0.01791
## 
## 
## $Edad_inicio
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.86818, p-value = 1.32e-09
## 
## 
## $Educacion
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.92555, p-value = 1.545e-06
## 
## 
## $Num_depresiones
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.89711, p-value = 3.444e-08
## 
## 
## $Intentossuicidio
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.83309, p-value = 4.445e-11
## 
## 
## $Internaciones
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.73011, p-value = 1.893e-14
## 
## 
## $YMRS
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.91715, p-value = 4.654e-07
## 
## 
## $MADRS
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.98508, p-value = 0.149
## 
## 
## $TMTA
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.9816, p-value = 0.0651
## 
## 
## $FE
## 
##  Shapiro-Wilk normality test
## 
## data:  lm(x ~ melanco$Melanco_SMPI)$residuals
## W = 0.97989, p-value = 0.04325

Metodo grafico

Tanto el metodo analitico como el grafico indican no cumplimiento del supuesto de normalidad para: Edad, Edad_inicio, Educación, Num_depresiones, Intentossuicidio, Iternaciones, YMRS y FE. Se cumple el supuesto solo para MADRS y TMTA.

Homocedasticidad

Metodo analitico
Test de Levene

lapply(melanco[,c("Edad", "Edad_inicio", "Educacion", "Num_depresiones",
                  "Intentossuicidio", "Internaciones", "YMRS", "MADRS",
                  "TMTA", "FE")],
       function(x) leveneTest(x~ melanco$Melanco_SMPI))

## $Edad
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   1  5.0989 0.02557 *
##       133                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Edad_inicio
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1   4e-04 0.9841
##       133               
## 
## $Educacion
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.5824 0.4467
##       133               
## 
## $Num_depresiones
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.0035 0.9528
##       133               
## 
## $Intentossuicidio
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group   1   7.249 0.008007 **
##       133                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Internaciones
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1   0.015 0.9029
##       133               
## 
## $YMRS
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)   
## group   1  9.0515 0.00314 **
##       133                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $MADRS
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  1.2917 0.2578
##       133               
## 
## $TMTA
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   1  4.1932 0.04256 *
##       133                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $FE
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.4498 0.5036
##       133

Se observa incumplimiento del supuesto de homocedasticidad para Edad, Intentossuicidio, YMRS y TMTA. Sin embargo, como se puede observar en la Figura 2, las distribuciones de estas variables entre los grupos difieren de una forma importante solo para Intentossuicidio e YMRS. Tanto Edad como TMTA tienen distribuciones similares en el metodo gráfico (Figura 2). Se cumple el supuesto para Edad_inicio, Educacion, Num_depresiones, Internaciones, MADRS y FE, lo que se condice con el analisis grafico de la Figura 2.

4.3.3 Resumen de supuestos y tabla 1

En resumen, luego de haber aplicado metodos analiticos y graficos, se decide realizar las comparaciones bivariadas que se muestran en la tabla 1 con las siguientes pruebas:

V. categóricas	Esperados >5	Prueba tabla 1
Psi_actuales	No	Test de Fisher
IFD_B_BDZ	No	Test de Fisher
IFD_B_AD	No	Test de Fiser
IFD_B_EA	No	Test de Fisher
IFD_B_AP	No	Test de Fisher
IFD_B_Ach	No	Test de Fisher
Sexo	Sí	Chi2
TBP	Sí	Chi2

V. continuas	Normalidad	Homocedasticidad	Prueba tabla 1
Edad	No	Sí	Wilcoxon
Edad_inicio	No	Sí	Wilcoxon
Educacion	No	Sí	Wilcoxon
Num_depresiones	No	Sí	Wilcoxon
Intentossuicidio	No	No	Test de la mediana
Internaciones	No	Sí	Wilcoxon
YMRS	No	No	Test de la mediana
MADRS	Sí	Sí	Prueba t
TMTA	Sí	Sí	Prueba t
FE	No	Sí	Wilcoxon

Tabla 1

	No melancolico (N=100)	Melancolico (N=35)	Total (N=135)	p value
Sexo				0.459
Hombre	25 (25.0%)	11 (31.4%)	36 (26.7%)
Mujer	75 (75.0%)	24 (68.6%)	99 (73.3%)
Edad				0.031
Mean	35.370	42.143	37.126
SD	12.368	15.677	13.574
Median	35.000	41.000	35.000
Q1, Q3	25.750, 43.250	29.500, 54.500	26.000, 47.000
Educacion				0.116
Mean	11.490	12.000	11.622
SD	2.368	3.039	2.556
Median	12.000	12.000	12.000
Q1, Q3	10.000, 12.000	12.000, 13.500	10.500, 12.000
TBP				0.103
TDM	56 (56.0%)	14 (40.0%)	70 (51.9%)
TBP	44 (44.0%)	21 (60.0%)	65 (48.1%)
Edad_inicio				0.004
Mean	22.130	27.171	23.437
SD	12.214	11.925	12.297
Median	18.000	24.000	20.000
Q1, Q3	13.000, 26.000	19.000, 31.000	14.000, 29.500
Psi_actuales				< 0.001
No	96 (96.0%)	26 (74.3%)	122 (90.4%)
Si	4 (4.0%)	9 (25.7%)	13 (9.6%)
Num_depresiones				0.121
Mean	5.300	4.486	5.089
SD	3.298	3.697	3.411
Median	5.000	3.000	4.000
Q1, Q3	2.000, 9.000	2.000, 7.500	2.000, 9.000
Internaciones				0.925
Mean	2.180	2.229	2.193
SD	2.451	2.340	2.414
Median	1.000	1.000	1.000
Q1, Q3	1.000, 3.000	1.000, 3.000	1.000, 3.000
Intentossuicidio				0.002
Mean	2.930	1.171	2.474
SD	2.996	1.618	2.809
Median	2.000	1.000	2.000
Q1, Q3	1.000, 4.000	0.000, 1.000	1.000, 3.000
MADRS				< 0.001
Mean	29.150	34.029	30.415
SD	6.908	8.194	7.543
Median	29.000	33.000	30.000
Q1, Q3	24.000, 33.250	29.000, 39.500	24.500, 36.000
YMRS				< 0.001
Mean	2.390	0.771	1.970
SD	1.989	1.374	1.977
Median	2.000	0.000	2.000
Q1, Q3	0.750, 4.000	0.000, 1.000	0.000, 4.000
IFD_B_BZD				0.009
0	10 (10.0%)	7 (20.0%)	17 (12.6%)
2	6 (6.0%)	5 (14.3%)	11 (8.1%)
3	7 (7.0%)	3 (8.6%)	10 (7.4%)
4	46 (46.0%)	18 (51.4%)	64 (47.4%)
5	31 (31.0%)	2 (5.7%)	33 (24.4%)
IFD_B_AD				0.582
0	48 (48.0%)	22 (62.9%)	70 (51.9%)
2	16 (16.0%)	3 (8.6%)	19 (14.1%)
3	17 (17.0%)	5 (14.3%)	22 (16.3%)
4	17 (17.0%)	4 (11.4%)	21 (15.6%)
5	2 (2.0%)	1 (2.9%)	3 (2.2%)
IFD_B_EA				0.403
0	57 (57.0%)	17 (48.6%)	74 (54.8%)
2	30 (30.0%)	9 (25.7%)	39 (28.9%)
3	7 (7.0%)	4 (11.4%)	11 (8.1%)
4	4 (4.0%)	3 (8.6%)	7 (5.2%)
5	2 (2.0%)	2 (5.7%)	4 (3.0%)
IFD_B_AP				0.146
0	12 (12.0%)	2 (5.7%)	14 (10.4%)
2	60 (60.0%)	21 (60.0%)	81 (60.0%)
3	18 (18.0%)	4 (11.4%)	22 (16.3%)
4	9 (9.0%)	5 (14.3%)	14 (10.4%)
5	1 (1.0%)	3 (8.6%)	4 (3.0%)
IFD_B_Ach				1.000
0	98 (98.0%)	34 (97.1%)	132 (97.8%)
2	2 (2.0%)	1 (2.9%)	3 (2.2%)
TMTA				0.002
Mean	0.273	-0.403	0.098
SD	0.994	1.376	1.140
Median	0.470	-0.260	0.370
Q1, Q3	-0.430, 0.935	-1.125, 0.520	-0.610, 0.835
FE				< 0.001
Mean	-0.305	-1.154	-0.525
SD	0.909	1.112	1.032
Median	-0.333	-1.000	-0.503
Q1, Q3	-0.781, 0.208	-1.447, -0.598	-1.020, 0.128

5 Datos faltantes

Cantidad de celdas con datos faltantes y casos perdidos por tener datos faltantes

sum(is.na(melanco))

## [1] 0

which(!complete.cases(melanco))

## integer(0)

Visualizacion de datos faltantes

vis_miss(melanco, sort=TRUE)

vis_miss(melanco, cluster = TRUE)

vis_dat(melanco)

En este caso, no tenemos datos faltantes.

6 Correlacion entre variables continuas

Cramos una nueva base que incluye solamente las variables continuas y luego realizamos las correlaciones.

variables_continuas<- melanco[,c("Edad", "Edad_inicio","Educacion", "Num_depresiones", "Intentossuicidio","Internaciones","YMRS", "MADRS", "TMTA", "FE")]

round(cor(variables_continuas),2)

##                   Edad Edad_inicio Educacion Num_depresiones Intentossuicidio
## Edad              1.00        0.60      0.10            0.02            -0.19
## Edad_inicio       0.60        1.00      0.25           -0.47            -0.29
## Educacion         0.10        0.25      1.00           -0.15            -0.16
## Num_depresiones   0.02       -0.47     -0.15            1.00             0.46
## Intentossuicidio -0.19       -0.29     -0.16            0.46             1.00
## Internaciones     0.10       -0.14      0.02            0.41             0.48
## YMRS             -0.19       -0.23     -0.14            0.06             0.16
## MADRS             0.06        0.02     -0.15            0.07            -0.03
## TMTA             -0.40       -0.20     -0.15           -0.03             0.07
## FE               -0.24       -0.14     -0.02           -0.05            -0.03
##                  Internaciones  YMRS MADRS  TMTA    FE
## Edad                      0.10 -0.19  0.06 -0.40 -0.24
## Edad_inicio              -0.14 -0.23  0.02 -0.20 -0.14
## Educacion                 0.02 -0.14 -0.15 -0.15 -0.02
## Num_depresiones           0.41  0.06  0.07 -0.03 -0.05
## Intentossuicidio          0.48  0.16 -0.03  0.07 -0.03
## Internaciones             1.00  0.09  0.10 -0.19 -0.17
## YMRS                      0.09  1.00  0.01  0.05  0.05
## MADRS                     0.10  0.01  1.00 -0.25 -0.26
## TMTA                     -0.19  0.05 -0.25  1.00  0.66
## FE                       -0.17  0.05 -0.26  0.66  1.00

corrplot(cor(variables_continuas), method = 'ellipse', order = 'AOE', type = 'upper')

Respecto de nuestra variable de interés FE, la misma tuvo correlaciones negativas con Edad y el puntaje en MADRS. Se observó Correlacion positiva de FE y TMTA.

Por otro lado, observamos correlacion positiva entre:
numero de depresiones con nro de intentos suicidas e internaciones;
intentos suicidas con internaciones;
edad y edad de inicio;
edad de inicio y educacion.

Correlacion negativa entre:
edad con intentos suicidas, puntaje YMRS y TMTA
edad de inicio con numero de depresiones, intentos suicidas y puntaje YMRS;
TMTA e internaciones; MADRS y TMTA.

7 Modelos lineales

Se proponen modelos para la variable respuesta FE. Se espera que una de las variables clínicas explicativas sea pertenecer al grupo de pacientes con melancolia (Melanco_SMPI = 1).

7.1 Modelo lineales simples

7.1.1 Modelos simples: Variables Categóricas Dicotómicas

	FE			FE			FE			FE
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.44	-0.78 – -0.10	0.011	-0.31	-0.55 – -0.07	0.012	-0.48	-0.67 – -0.30	<0.001	-0.30	-0.50 – -0.11	0.002
Sexo [Mujer]	-0.11	-0.51 – 0.29	0.583
TBP [TBP]				-0.45	-0.79 – -0.11	0.011
Psi actuales [Si]							-0.44	-1.03 – 0.15	0.143
Melanco SMPI [Melancolico]										-0.85	-1.22 – -0.47	<0.001
Observations	135			135			135			135
R² / R² adjusted	0.002 / -0.005			0.048 / 0.041			0.016 / 0.009			0.131 / 0.125

7.1.2 Modelos simples: Variables Categóricas Ordinales

## $emmeans
##  IFD_B_BZD emmean    SE  df lower.CL upper.CL
##  0         -0.485 0.249 130   -0.979  0.00823
##  2         -1.114 0.310 130   -1.727 -0.50040
##  3         -0.724 0.325 130   -1.368 -0.08085
##  4         -0.408 0.129 130   -0.662 -0.15377
##  5         -0.515 0.179 130   -0.869 -0.16073
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate    SE  df t.ratio p.value
##  IFD_B_BZD2 - IFD_B_BZD0  -0.6286 0.398 130  -1.580  0.3304
##  IFD_B_BZD3 - IFD_B_BZD0  -0.2390 0.410 130  -0.583  0.9020
##  IFD_B_BZD4 - IFD_B_BZD0   0.0772 0.281 130   0.275  0.9833
##  IFD_B_BZD5 - IFD_B_BZD0  -0.0297 0.307 130  -0.097  0.9986
## 
## P value adjustment: dunnettx method for 4 tests

## $emmeans
##  IFD_B_BZD emmean    SE  df lower.CL upper.CL
##  0         -0.485 0.249 130   -0.979  0.00823
##  2         -1.114 0.310 130   -1.727 -0.50040
##  3         -0.724 0.325 130   -1.368 -0.08085
##  4         -0.408 0.129 130   -0.662 -0.15377
##  5         -0.515 0.179 130   -0.869 -0.16073
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                estimate    SE  df lower.CL upper.CL
##  IFD_B_BZD2 - IFD_B_BZD0  -0.6286 0.398 130   -1.617    0.360
##  IFD_B_BZD3 - IFD_B_BZD0  -0.2390 0.410 130   -1.257    0.779
##  IFD_B_BZD4 - IFD_B_BZD0   0.0772 0.281 130   -0.620    0.774
##  IFD_B_BZD5 - IFD_B_BZD0  -0.0297 0.307 130   -0.792    0.733
## 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 4 estimates

## $emmeans
##  IFD_B_AD emmean    SE  df lower.CL upper.CL
##  0        -0.578 0.124 130   -0.823  -0.3321
##  2        -0.319 0.238 130   -0.790   0.1521
##  3        -0.357 0.221 130   -0.795   0.0804
##  4        -0.734 0.226 130   -1.182  -0.2856
##  5        -0.368 0.599 130   -1.553   0.8177
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_AD2 - IFD_B_AD0    0.259 0.268 130   0.963  0.7075
##  IFD_B_AD3 - IFD_B_AD0    0.220 0.254 130   0.868  0.7640
##  IFD_B_AD4 - IFD_B_AD0   -0.156 0.258 130  -0.605  0.8936
##  IFD_B_AD5 - IFD_B_AD0    0.210 0.612 130   0.343  0.9718
## 
## P value adjustment: dunnettx method for 4 tests

## $emmeans
##  IFD_B_AD emmean    SE  df lower.CL upper.CL
##  0        -0.578 0.124 130   -0.823  -0.3321
##  2        -0.319 0.238 130   -0.790   0.1521
##  3        -0.357 0.221 130   -0.795   0.0804
##  4        -0.734 0.226 130   -1.182  -0.2856
##  5        -0.368 0.599 130   -1.553   0.8177
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df lower.CL upper.CL
##  IFD_B_AD2 - IFD_B_AD0    0.259 0.268 130   -0.408    0.926
##  IFD_B_AD3 - IFD_B_AD0    0.220 0.254 130   -0.410    0.850
##  IFD_B_AD4 - IFD_B_AD0   -0.156 0.258 130   -0.798    0.485
##  IFD_B_AD5 - IFD_B_AD0    0.210 0.612 130   -1.310    1.730
## 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 4 estimates

## $emmeans
##  IFD_B_AP emmean    SE  df lower.CL upper.CL
##  0        -0.168 0.274 130   -0.710    0.375
##  2        -0.543 0.114 130   -0.768   -0.317
##  3        -0.606 0.219 130   -1.039   -0.173
##  4        -0.381 0.274 130   -0.924    0.161
##  5        -1.468 0.513 130   -2.483   -0.453
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_AP2 - IFD_B_AP0   -0.375 0.297 130  -1.263  0.5170
##  IFD_B_AP3 - IFD_B_AP0   -0.439 0.351 130  -1.250  0.5253
##  IFD_B_AP4 - IFD_B_AP0   -0.214 0.388 130  -0.551  0.9141
##  IFD_B_AP5 - IFD_B_AP0   -1.301 0.582 130  -2.236  0.0911
## 
## P value adjustment: dunnettx method for 4 tests

## $emmeans
##  IFD_B_AP emmean    SE  df lower.CL upper.CL
##  0        -0.168 0.274 130   -0.710    0.375
##  2        -0.543 0.114 130   -0.768   -0.317
##  3        -0.606 0.219 130   -1.039   -0.173
##  4        -0.381 0.274 130   -0.924    0.161
##  5        -1.468 0.513 130   -2.483   -0.453
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df lower.CL upper.CL
##  IFD_B_AP2 - IFD_B_AP0   -0.375 0.297 130    -1.11    0.363
##  IFD_B_AP3 - IFD_B_AP0   -0.439 0.351 130    -1.31    0.433
##  IFD_B_AP4 - IFD_B_AP0   -0.214 0.388 130    -1.18    0.750
##  IFD_B_AP5 - IFD_B_AP0   -1.301 0.582 130    -2.75    0.144
## 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 4 estimates

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.287 0.115 130   -0.514  -0.0592
##  2        -0.717 0.158 130   -1.031  -0.4039
##  3        -0.959 0.298 130   -1.549  -0.3684
##  4        -0.500 0.374 130   -1.240   0.2400
##  5        -1.903 0.495 130   -2.882  -0.9244
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_EA2 - IFD_B_EA0   -0.431 0.196 130  -2.199  0.0991
##  IFD_B_EA3 - IFD_B_EA0   -0.672 0.320 130  -2.101  0.1232
##  IFD_B_EA4 - IFD_B_EA0   -0.213 0.391 130  -0.545  0.9162
##  IFD_B_EA5 - IFD_B_EA0   -1.617 0.508 130  -3.182  0.0069
## 
## P value adjustment: dunnettx method for 4 tests

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.287 0.115 130   -0.514  -0.0592
##  2        -0.717 0.158 130   -1.031  -0.4039
##  3        -0.959 0.298 130   -1.549  -0.3684
##  4        -0.500 0.374 130   -1.240   0.2400
##  5        -1.903 0.495 130   -2.882  -0.9244
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df lower.CL upper.CL
##  IFD_B_EA2 - IFD_B_EA0   -0.431 0.196 130   -0.917   0.0558
##  IFD_B_EA3 - IFD_B_EA0   -0.672 0.320 130   -1.467   0.1224
##  IFD_B_EA4 - IFD_B_EA0   -0.213 0.391 130   -1.185   0.7590
##  IFD_B_EA5 - IFD_B_EA0   -1.617 0.508 130   -2.879  -0.3545
## 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 4 estimates

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.287 0.115 130   -0.514  -0.0592
##  2        -0.717 0.158 130   -1.031  -0.4039
##  3        -0.959 0.298 130   -1.549  -0.3684
##  4        -0.500 0.374 130   -1.240   0.2400
##  5        -1.903 0.495 130   -2.882  -0.9244
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_EA0 - IFD_B_EA2    0.431 0.196 130   2.199  0.1865
##  IFD_B_EA0 - IFD_B_EA3    0.672 0.320 130   2.101  0.2258
##  IFD_B_EA0 - IFD_B_EA4    0.213 0.391 130   0.545  0.9824
##  IFD_B_EA0 - IFD_B_EA5    1.617 0.508 130   3.182  0.0155
##  IFD_B_EA2 - IFD_B_EA3    0.241 0.338 130   0.714  0.9529
##  IFD_B_EA2 - IFD_B_EA4   -0.217 0.406 130  -0.535  0.9835
##  IFD_B_EA2 - IFD_B_EA5    1.186 0.520 130   2.282  0.1573
##  IFD_B_EA3 - IFD_B_EA4   -0.459 0.478 130  -0.959  0.8729
##  IFD_B_EA3 - IFD_B_EA5    0.945 0.578 130   1.635  0.4782
##  IFD_B_EA4 - IFD_B_EA5    1.403 0.620 130   2.262  0.1640
## 
## P value adjustment: tukey method for comparing a family of 5 estimates

## 
## Call:
## lm(formula = FE ~ IFD_B_Ach, data = melanco)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.13773 -0.50773  0.05561  0.64061  2.62227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.51227    0.08984  -5.702 7.31e-08 ***
## IFD_B_Ach2  -0.56773    0.60267  -0.942    0.348    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.032 on 133 degrees of freedom
## Multiple R-squared:  0.006628,   Adjusted R-squared:  -0.000841 
## F-statistic: 0.8874 on 1 and 133 DF,  p-value: 0.3479

	FE			FE			FE			FE			FE
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.49	-0.98 – 0.01	0.054	-0.58	-0.82 – -0.33	<0.001	-0.17	-0.71 – 0.37	0.542	-0.29	-0.51 – -0.06	0.014	-0.51	-0.69 – -0.33	<0.001
IFD B BZD [2]	-0.63	-1.42 – 0.16	0.117
IFD B BZD [3]	-0.24	-1.05 – 0.57	0.561
IFD B BZD [4]	0.08	-0.48 – 0.63	0.784
IFD B BZD [5]	-0.03	-0.64 – 0.58	0.923
IFD B AD [2]				0.26	-0.27 – 0.79	0.337
IFD B AD [3]				0.22	-0.28 – 0.72	0.387
IFD B AD [4]				-0.16	-0.67 – 0.35	0.547
IFD B AD [5]				0.21	-1.00 – 1.42	0.732
IFD B AP [2]							-0.38	-0.96 – 0.21	0.209
IFD B AP [3]							-0.44	-1.13 – 0.26	0.213
IFD B AP [4]							-0.21	-0.98 – 0.55	0.583
IFD B AP [5]							-1.30	-2.45 – -0.15	0.027
IFD B EA [2]										-0.43	-0.82 – -0.04	0.030
IFD B EA [3]										-0.67	-1.30 – -0.04	0.038
IFD B EA [4]										-0.21	-0.99 – 0.56	0.587
IFD B EA [5]										-1.62	-2.62 – -0.61	0.002
IFD B Ach [2]													-0.57	-1.76 – 0.62	0.348
Observations	135			135			135			135			135
R² / R² adjusted	0.036 / 0.006			0.018 / -0.012			0.041 / 0.011			0.107 / 0.080			0.007 / -0.001

7.1.3 Modelos simples: Variables continuas

	FE			FE			FE			FE			FE			FE			FE
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.25	-0.63 – 0.13	0.193	-0.45	-0.77 – -0.14	0.005	-0.37	-0.60 – -0.13	0.002	-0.50	-0.73 – -0.26	<0.001	0.56	-0.15 – 1.27	0.122	-0.58	-0.83 – -0.33	<0.001	-0.58	-0.72 – -0.45	<0.001
Edad inicio	-0.01	-0.03 – 0.00	0.105
Num depresiones				-0.01	-0.07 – 0.04	0.592
Internaciones							-0.07	-0.14 – 0.00	0.051
Intentossuicidio										-0.01	-0.07 – 0.05	0.747
MADRS													-0.04	-0.06 – -0.01	0.002
YMRS																0.03	-0.06 – 0.12	0.529
TMTA																			0.59	0.48 – 0.71	<0.001
Observations	135			135			135			135			135			135			135
R² / R² adjusted	0.020 / 0.012			0.002 / -0.005			0.028 / 0.021			0.001 / -0.007			0.068 / 0.061			0.003 / -0.005			0.430 / 0.425

Los modelos lineales simples sugieren que de las 19 variables, solo fueron significativas: TBP, Melanco_SMPI, IFD_B_EA, MADRS, y TMTA. Los dos modelos simples con R2 mas altos fueron Melanco_SMPI (R2 = 0.13) y TMTA (R2 = 0.43).

De estas variables, las que tambien habian mostrado significacion en el modelo con todas las variables fueron: Melanco_SMPI, IFD_B_EA, y TMTA. (Yo esto lo plantearia abajo al reves, como cambian cuando las pones en un modelo multiple)

7.2 Modelo con todas las variables

Comenzamos proponiendo el modelo que contiene a todas las variables (modelo1), con excepción de las variables edad y educación, ya que la variable respuesta se encuentra estandarizada por esas variables. Además, calculamos los VIF y comparamos los grupos en las variables con más de una categoría con el comando emmeans

melanco2<-melanco[-c(2,3)]
modelo1 <-lm(FE ~., data = melanco2)
summary(modelo1)

## 
## Call:
## lm(formula = FE ~ ., data = melanco2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.65318 -0.42340 -0.02862  0.40753  1.74254 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.777088   0.494723   1.571  0.11922    
## SexoMujer               -0.210169   0.171714  -1.224  0.22369    
## TBPTBP                  -0.150792   0.153942  -0.980  0.32955    
## Edad_inicio             -0.007842   0.006913  -1.134  0.25918    
## Num_depresiones          0.005910   0.026588   0.222  0.82452    
## Melanco_SMPIMelancolico -0.463897   0.201802  -2.299  0.02348 *  
## Intentossuicidio        -0.069111   0.030414  -2.272  0.02508 *  
## Internaciones            0.046656   0.036031   1.295  0.19817    
## Psi_actualesSi           0.452512   0.276313   1.638  0.10445    
## YMRS                     0.012366   0.039235   0.315  0.75325    
## MADRS                   -0.014699   0.010049  -1.463  0.14650    
## IFD_B_BZD2               0.214398   0.318836   0.672  0.50277    
## IFD_B_BZD3               0.064165   0.339418   0.189  0.85042    
## IFD_B_BZD4               0.029545   0.230380   0.128  0.89820    
## IFD_B_BZD5              -0.227920   0.266974  -0.854  0.39519    
## IFD_B_AD2                0.102444   0.220501   0.465  0.64317    
## IFD_B_AD3                0.282108   0.222946   1.265  0.20852    
## IFD_B_AD4                0.060013   0.227298   0.264  0.79227    
## IFD_B_AD5               -0.332880   0.474070  -0.702  0.48411    
## IFD_B_EA2               -0.213481   0.167448  -1.275  0.20513    
## IFD_B_EA3               -0.784214   0.275747  -2.844  0.00535 ** 
## IFD_B_EA4               -0.230971   0.326574  -0.707  0.48096    
## IFD_B_EA5               -0.651822   0.419980  -1.552  0.12364    
## IFD_B_AP2               -0.294038   0.250238  -1.175  0.24261    
## IFD_B_AP3               -0.378033   0.298542  -1.266  0.20819    
## IFD_B_AP4               -0.173286   0.321356  -0.539  0.59086    
## IFD_B_AP5               -1.218439   0.525209  -2.320  0.02226 *  
## IFD_B_Ach2              -0.102882   0.493732  -0.208  0.83534    
## TMTA                     0.562325   0.070495   7.977 1.87e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7247 on 106 degrees of freedom
## Multiple R-squared:  0.6098, Adjusted R-squared:  0.5067 
## F-statistic: 5.915 on 28 and 106 DF,  p-value: 8.513e-12

Anova(modelo1, type="III")

vif_values_FE<-vif(modelo1)[,1]
vif_values_FE

##             Sexo              TBP      Edad_inicio  Num_depresiones 
##         1.482256         1.520897         1.843754         2.098524 
##     Melanco_SMPI Intentossuicidio    Internaciones     Psi_actuales 
##         2.010431         1.862943         1.930672         1.707963 
##             YMRS            MADRS        IFD_B_BZD         IFD_B_AD 
##         1.535594         1.466074         2.993957         3.250313 
##         IFD_B_EA         IFD_B_AP        IFD_B_Ach             TMTA 
##         2.870401         3.813509         1.361599         1.647061

barplot(vif_values_FE, axes = T, main = "VIF Values", col = "steelblue", las=2)

comp_ea_modelo1<-emmeans(modelo1, specs = trt.vs.ctrl ~ IFD_B_EA, ref = 1, rg.limit =40000)
comp_ea_modelo1 #continua siendo significativo nivel 0 vs 3 de ea

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.425 0.299 106    -1.02   0.1679
##  2        -0.639 0.349 106    -1.33   0.0534
##  3        -1.209 0.352 106    -1.91  -0.5124
##  4        -0.656 0.410 106    -1.47   0.1566
##  5        -1.077 0.505 106    -2.08  -0.0765
## 
## Results are averaged over the levels of: Sexo, TBP, Melanco_SMPI, Psi_actuales, IFD_B_BZD, IFD_B_AD, IFD_B_AP, IFD_B_Ach 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_EA2 - IFD_B_EA0   -0.213 0.167 106  -1.275  0.5104
##  IFD_B_EA3 - IFD_B_EA0   -0.784 0.276 106  -2.844  0.0195
##  IFD_B_EA4 - IFD_B_EA0   -0.231 0.327 106  -0.707  0.8485
##  IFD_B_EA5 - IFD_B_EA0   -0.652 0.420 106  -1.552  0.3461
## 
## Results are averaged over the levels of: Sexo, TBP, Melanco_SMPI, Psi_actuales, IFD_B_BZD, IFD_B_AD, IFD_B_AP, IFD_B_Ach 
## P value adjustment: dunnettx method for 4 tests

En el modelo en el que busco explicar el desempeño en FE a partir de todas las covariables, se observan signficativas mi variable de interes, Melanco_SMPI y resultan significativas también Intentos de suicidio y TMTA, IFD_B_EA da un p-valor limitrofe (0.055). El modelo obtiene un R2 ajustado de .507. No se observan valores altos de VIF para ninguna variable. Es un modelo con 29 parámetros.

Supuestos

Supuesto de normalidad

qqnorm (modelo1$residuals)
qqline(modelo1$residuals)

shapiro.test (modelo1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo1$residuals
## W = 0.99055, p-value = 0.4952

No rechazo la normalidad de los residuos.

Supuesto de homocedasticidad

plot(modelo1$fitted.values, modelo1$residuals)
abline(0,0)

No rechazo el supuesto de igualdad de varianzas.

7.3 Modelos por seleccion automatica de variables

Probamos con metodos de seleccion automatica para ver la coincidencia entre el modelo con todas las variables, los modelos simples y la importancia clinica.

mat_modelo<-model.matrix(modelo1) 
Xy<-data.frame(cbind(mat_modelo[,-1], melanco2$FE))

reg.bestglm1<-bestglm (Xy=Xy, 
                      IC="AIC", 
                      method= "forward")
# reg.bestglm1$Subsets

reg.bestglm2<-bestglm (Xy=Xy, 
                      IC="AIC", 
                      method= "backward")
# reg.bestglm2$Subsets

reg.bestglm3<-bestglm (Xy=Xy, 
                      IC="AIC", 
                      method= "exhaustive")
# reg.bestglm3$Subsets

reg.bestglm1$BestModel

## 
## Call:
## lm(formula = y ~ ., data = data.frame(Xy[, c(bestset[-1], FALSE), 
##     drop = FALSE], y = y))
## 
## Coefficients:
##             (Intercept)                   TBPTBP  Melanco_SMPIMelancolico  
##                 0.25917                 -0.20443                 -0.48078  
##        Intentossuicidio           Psi_actualesSi                    MADRS  
##                -0.03806                  0.44807                 -0.01532  
##              IFD_B_BZD5                IFD_B_AD3                IFD_B_EA3  
##                -0.23677                  0.28240                 -0.57565  
##               IFD_B_EA5                IFD_B_AP5                     TMTA  
##                -0.51615                 -0.80125                  0.55141

reg.bestglm2$BestModel

## 
## Call:
## lm(formula = y ~ ., data = data.frame(Xy[, c(bestset[-1], FALSE), 
##     drop = FALSE], y = y))
## 
## Coefficients:
##             (Intercept)                   TBPTBP  Melanco_SMPIMelancolico  
##                 0.25917                 -0.20443                 -0.48078  
##        Intentossuicidio           Psi_actualesSi                    MADRS  
##                -0.03806                  0.44807                 -0.01532  
##              IFD_B_BZD5                IFD_B_AD3                IFD_B_EA3  
##                -0.23677                  0.28240                 -0.57565  
##               IFD_B_EA5                IFD_B_AP5                     TMTA  
##                -0.51615                 -0.80125                  0.55141

reg.bestglm3$BestModel

## 
## Call:
## lm(formula = y ~ ., data = data.frame(Xy[, c(bestset[-1], FALSE), 
##     drop = FALSE], y = y))
## 
## Coefficients:
##             (Intercept)                   TBPTBP  Melanco_SMPIMelancolico  
##                 0.25917                 -0.20443                 -0.48078  
##        Intentossuicidio           Psi_actualesSi                    MADRS  
##                -0.03806                  0.44807                 -0.01532  
##              IFD_B_BZD5                IFD_B_AD3                IFD_B_EA3  
##                -0.23677                  0.28240                 -0.57565  
##               IFD_B_EA5                IFD_B_AP5                     TMTA  
##                -0.51615                 -0.80125                  0.55141

Todos los metodos stepwise seleccionan las mismas variables: TBP, Melanco_SMPI, Intentossuicidio,Psi_actuales, MADRS,IFD_B_BZD, IFD_B_AD, IFD_B_EA, IFD_B_AP, TMTA.

7.3.1 Random Forest

set.seed(44)
rf<-randomForest(FE~., data = melanco2, importance = TRUE)


rf$importance

##                        %IncMSE IncNodePurity
## Sexo             -0.0082070911      2.205236
## TBP               0.0068838230      2.112003
## Edad_inicio       0.0478652017     14.119618
## Num_depresiones   0.0280036185      5.499344
## Melanco_SMPI      0.0746667911      7.825772
## Intentossuicidio  0.0062102283      4.447652
## Internaciones     0.0034520396      4.718148
## Psi_actuales     -0.0009904466      1.232982
## YMRS              0.0097158226      3.761420
## MADRS             0.0527219393     12.738989
## IFD_B_BZD        -0.0067685963      7.782563
## IFD_B_AD         -0.0016124690      5.070561
## IFD_B_EA          0.0506893093     10.690487
## IFD_B_AP         -0.0046515477      5.521846
## IFD_B_Ach        -0.0015934187      0.214509
## TMTA              0.4523593629     44.640935

varImpPlot(rf) #Aca puedo ver cuanto aumenta el MSE si saco la variable, parece que dsp de IFD_EA no hay cambios importantes. Las variables seleccionadas por Random Forest serian TMTA, Melanco_SMPI, MADRS, IFD_B_EA

v<-as.vector(rf$importance[,1])
w<-(as.vector((colnames(melanco2[1:16]))))
DF<-cbind(w,v)
DF<-as.data.frame(DF)
str(DF)

## 'data.frame':    16 obs. of  2 variables:
##  $ w: chr  "Sexo" "TBP" "Edad_inicio" "Num_depresiones" ...
##  $ v: chr  "-0.00820709109253545" "0.00688382302393686" "0.0478652016563371" "0.028003618454361" ...

DF<-DF %>% mutate(v=as.numeric(v),
              w=as.factor(w))

ggplot(DF, aes(x=reorder(w,v), y=v,fill=w))+ 
  geom_bar(stat="identity", position="dodge")+ coord_flip()+
  ylab("Variable Importance")+
  xlab("")+
  theme(legend.position = "none")

7.3.1.1 Ajuste de modelo con las variables de selección automática

Ajustamos un modelo con las variables:TBP, Melanco_SMPI, Intentossuicidio,Psi_actuales, MADRS,IFD_B_BZD, IFD_B_AD, IFD_B_EA, IFD_B_AP, TMTA . Lo comparamos con el modelo que contenía todas las variables.

## 
## Call:
## lm(formula = FE ~ TBP + Melanco_SMPI + Intentossuicidio + Psi_actuales + 
##     MADRS + IFD_B_BZD + IFD_B_AD + IFD_B_EA + IFD_B_AP + TMTA, 
##     data = melanco2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.63353 -0.41561 -0.06449  0.44165  1.56702 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.373054   0.356522   1.046  0.29764    
## TBPTBP                  -0.169040   0.148272  -1.140  0.25669    
## Melanco_SMPIMelancolico -0.493489   0.176985  -2.788  0.00623 ** 
## Intentossuicidio        -0.040272   0.024550  -1.640  0.10373    
## Psi_actualesSi           0.439342   0.270002   1.627  0.10651    
## MADRS                   -0.013592   0.009655  -1.408  0.16197    
## IFD_B_BZD2               0.203086   0.311217   0.653  0.51538    
## IFD_B_BZD3               0.070128   0.320541   0.219  0.82722    
## IFD_B_BZD4               0.044785   0.223571   0.200  0.84160    
## IFD_B_BZD5              -0.162133   0.256033  -0.633  0.52786    
## IFD_B_AD2                0.015940   0.209168   0.076  0.93939    
## IFD_B_AD3                0.237365   0.212267   1.118  0.26586    
## IFD_B_AD4               -0.016657   0.215187  -0.077  0.93844    
## IFD_B_AD5               -0.329958   0.453359  -0.728  0.46825    
## IFD_B_EA2               -0.139327   0.159435  -0.874  0.38405    
## IFD_B_EA3               -0.616729   0.259733  -2.374  0.01928 *  
## IFD_B_EA4               -0.006956   0.304778  -0.023  0.98183    
## IFD_B_EA5               -0.603185   0.403127  -1.496  0.13740    
## IFD_B_AP2               -0.181666   0.238375  -0.762  0.44760    
## IFD_B_AP3               -0.317913   0.286426  -1.110  0.26941    
## IFD_B_AP4               -0.066312   0.311128  -0.213  0.83161    
## IFD_B_AP5               -1.022451   0.496225  -2.060  0.04167 *  
## TMTA                     0.557726   0.066371   8.403 1.54e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7209 on 112 degrees of freedom
## Multiple R-squared:  0.592,  Adjusted R-squared:  0.5118 
## F-statistic: 7.386 on 22 and 112 DF,  p-value: 1.999e-13

El modelo con las variables seleccionadas automaticamente tiene un R2 = 0.512. No mejora mucho respecto del que contenia todas las variables. Sigue siendo un modelo con muchos parámetros (23). Son significativas Melanco_SMPI y TMTA

8 Modelo final unificando resultados y usando criterios clínicos

Justificacion

Consideramos las variables que de forma repetida mostraron relevancia algunos de los modelos analizados. Estas son: Melanco_SMPI (nuestra variable de mayor interes), Intentossuicidio, IFD_B_EA, y TMTA.

Además, sumamos otras variables que tienen importancia clinica, que son relevantes para el desempeño de las otras variables y/o que mostraron diferencias en la tabla 1. Además estas variables fueron seleccionadas en los métodos automáticos de selección de variables: MADRS, Psi_actuales y TBP. Nos parece importante que los modelos contengan a la variable MADRS (medida de severidad de la depresion) y sintomas psicoticos actuales (otra medida de severidad de la depresion). Depresiones más severas se asocian clínicamente a peor desempeño en funciones ejecutivas. Queremos evaluar si controlando por severidad de la depresión, el diagnóstico de melancolía continúa siendo un predictor de peor funcionamiento ejecutivo

8.1 Planteo del modelo final

El modelo final incluye: Melanco_SMPI + Intentossuicidio + Psi_actuales + MADRS + TBP + IFD_B_EA + TMTA

modelo_final <- lm(FE ~ Melanco_SMPI + Intentossuicidio + Psi_actuales +
                     MADRS + TBP + IFD_B_EA  + TMTA,
                   data = melanco2)
summary(modelo_final)

## 
## Call:
## lm(formula = FE ~ Melanco_SMPI + Intentossuicidio + Psi_actuales + 
##     MADRS + TBP + IFD_B_EA + TMTA, data = melanco2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.42347 -0.43769 -0.04789  0.49984  1.67192 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.213273   0.288992   0.738  0.46191    
## Melanco_SMPIMelancolico -0.486630   0.166103  -2.930  0.00404 ** 
## Intentossuicidio        -0.045451   0.023350  -1.947  0.05385 .  
## Psi_actualesSi           0.237212   0.227131   1.044  0.29834    
## MADRS                   -0.011520   0.009094  -1.267  0.20763    
## TBPTBP                  -0.227897   0.135209  -1.686  0.09440 .  
## IFD_B_EA2               -0.130144   0.151430  -0.859  0.39176    
## IFD_B_EA3               -0.669192   0.239628  -2.793  0.00606 ** 
## IFD_B_EA4                0.015105   0.297303   0.051  0.95956    
## IFD_B_EA5               -0.733784   0.381325  -1.924  0.05661 .  
## TMTA                     0.519832   0.058754   8.848 7.55e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7194 on 124 degrees of freedom
## Multiple R-squared:  0.5502, Adjusted R-squared:  0.5139 
## F-statistic: 15.17 on 10 and 124 DF,  p-value: < 2.2e-16

tab_model(modelo_final)

	FE
Predictors	Estimates	CI	p
(Intercept)	0.21	-0.36 – 0.79	0.462
Melanco SMPI [Melancolico]	-0.49	-0.82 – -0.16	0.004
Intentossuicidio	-0.05	-0.09 – 0.00	0.054
Psi actuales [Si]	0.24	-0.21 – 0.69	0.298
MADRS	-0.01	-0.03 – 0.01	0.208
TBP [TBP]	-0.23	-0.50 – 0.04	0.094
IFD B EA [2]	-0.13	-0.43 – 0.17	0.392
IFD B EA [3]	-0.67	-1.14 – -0.19	0.006
IFD B EA [4]	0.02	-0.57 – 0.60	0.960
IFD B EA [5]	-0.73	-1.49 – 0.02	0.057
TMTA	0.52	0.40 – 0.64	<0.001
Observations	135
R² / R² adjusted	0.550 / 0.514

Anova(modelo_final, type="III")

comp_ea_modelo_final<-emmeans(modelo_final, pairwise ~ IFD_B_EA, rg.limit =40000)
comp_ea_modelo_final #continua siendo significativo nivel 0 vs 3 de ea

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.437 0.123 124   -0.682   -0.193
##  2        -0.568 0.157 124   -0.878   -0.258
##  3        -1.107 0.230 124   -1.562   -0.652
##  4        -0.422 0.299 124   -1.013    0.169
##  5        -1.171 0.375 124   -1.912   -0.430
## 
## Results are averaged over the levels of: Melanco_SMPI, Psi_actuales, TBP 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_EA0 - IFD_B_EA2   0.1301 0.151 124   0.859  0.9110
##  IFD_B_EA0 - IFD_B_EA3   0.6692 0.240 124   2.793  0.0469
##  IFD_B_EA0 - IFD_B_EA4  -0.0151 0.297 124  -0.051  1.0000
##  IFD_B_EA0 - IFD_B_EA5   0.7338 0.381 124   1.924  0.3100
##  IFD_B_EA2 - IFD_B_EA3   0.5390 0.255 124   2.114  0.2208
##  IFD_B_EA2 - IFD_B_EA4  -0.1452 0.304 124  -0.478  0.9892
##  IFD_B_EA2 - IFD_B_EA5   0.6036 0.384 124   1.572  0.5179
##  IFD_B_EA3 - IFD_B_EA4  -0.6843 0.352 124  -1.942  0.3010
##  IFD_B_EA3 - IFD_B_EA5   0.0646 0.429 124   0.150  0.9999
##  IFD_B_EA4 - IFD_B_EA5   0.7489 0.461 124   1.626  0.4840
## 
## Results are averaged over the levels of: Melanco_SMPI, Psi_actuales, TBP 
## P value adjustment: tukey method for comparing a family of 5 estimates

comp_ea_2<-emmeans(modelo_final, specs = trt.vs.ctrl ~ IFD_B_EA, ref = 1)
comp_ea_2

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.437 0.123 124   -0.682   -0.193
##  2        -0.568 0.157 124   -0.878   -0.258
##  3        -1.107 0.230 124   -1.562   -0.652
##  4        -0.422 0.299 124   -1.013    0.169
##  5        -1.171 0.375 124   -1.912   -0.430
## 
## Results are averaged over the levels of: Melanco_SMPI, Psi_actuales, TBP 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  IFD_B_EA2 - IFD_B_EA0  -0.1301 0.151 124  -0.859  0.7687
##  IFD_B_EA3 - IFD_B_EA0  -0.6692 0.240 124  -2.793  0.0220
##  IFD_B_EA4 - IFD_B_EA0   0.0151 0.297 124   0.051  0.9997
##  IFD_B_EA5 - IFD_B_EA0  -0.7338 0.381 124  -1.924  0.1781
## 
## Results are averaged over the levels of: Melanco_SMPI, Psi_actuales, TBP 
## P value adjustment: dunnettx method for 4 tests

confint(comp_ea_2)

## $emmeans
##  IFD_B_EA emmean    SE  df lower.CL upper.CL
##  0        -0.437 0.123 124   -0.682   -0.193
##  2        -0.568 0.157 124   -0.878   -0.258
##  3        -1.107 0.230 124   -1.562   -0.652
##  4        -0.422 0.299 124   -1.013    0.169
##  5        -1.171 0.375 124   -1.912   -0.430
## 
## Results are averaged over the levels of: Melanco_SMPI, Psi_actuales, TBP 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df lower.CL upper.CL
##  IFD_B_EA2 - IFD_B_EA0  -0.1301 0.151 124   -0.507   0.2463
##  IFD_B_EA3 - IFD_B_EA0  -0.6692 0.240 124   -1.265  -0.0735
##  IFD_B_EA4 - IFD_B_EA0   0.0151 0.297 124   -0.724   0.7541
##  IFD_B_EA5 - IFD_B_EA0  -0.7338 0.381 124   -1.682   0.2141
## 
## Results are averaged over the levels of: Melanco_SMPI, Psi_actuales, TBP 
## Confidence level used: 0.95 
## Conf-level adjustment: dunnettx method for 4 estimates

Medidas de bondad de ajuste para el modelo final, aparentes y por k-fold Cross-Validation : MSE (mean squared error)- RMSE (root mean squared error)- MAE (mean absolute error)

ecm_modelo_final<-mean((predict(modelo_final, melanco2)- melanco2$FE)^2)
ecm_modelo_final

## [1] 0.4753049

rmse_modelo_final<-sqrt(ecm_modelo_final)
rmse_modelo_final

## [1] 0.6894236

mae_modelo_final<-mean(abs((predict(modelo_final, melanco2)- melanco2$FE)))
mae_modelo_final

## [1] 0.5330675

set.seed(42)
train_control <- trainControl(method = "cv",
                              number = 5)

ECM_modelo_final <- train(FE ~ Melanco_SMPI + Intentossuicidio + Psi_actuales +
                     MADRS + TBP + IFD_B_EA  + TMTA, data = melanco2,
               trControl = train_control,
               method = "lm")

print(ECM_modelo_final)

## Linear Regression 
## 
## 135 samples
##   7 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 107, 109, 108, 108, 108 
## Resampling results:
## 
##   RMSE       Rsquared  MAE      
##   0.7630343  0.485314  0.6045929
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

Supuestos
Normalidad y homocedasticidad.

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_final$residuals
## W = 0.98691, p-value = 0.2281

No tenemos evidencia para rechazar los supuestos de normalidad y homocedasticidad. El modelo final tiene un R2 ajustado de 0.51. En este modelo son significativas las variables Melanco_SMPI, IFD_B_EA, TMTA, y tiene un p valor limitrofe intentossucidio. Mientras que los coeficientes betas estimados de las variables IFD_B_EA, TMTA no se modifica respecto a los obtenidos en los modelos bivariados, el de Melanco_SMPI se reduce de -0.85 a -0.49.

Debido a que la variable TMT-A tiene una correlación fuerte con la variable respuesta (r=0.77). Ajustamos un modelo similar pero sin esta variable y obervamos el comportamiento de los coeficientes.

modelo_sin_TMTA <- lm(FE ~ Melanco_SMPI + Intentossuicidio + Psi_actuales +
                     MADRS + TBP + IFD_B_EA,
                   data = melanco2)
summary(modelo_sin_TMTA)

## 
## Call:
## lm(formula = FE ~ Melanco_SMPI + Intentossuicidio + Psi_actuales + 
##     MADRS + TBP + IFD_B_EA, data = melanco2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.52644 -0.39490 -0.03217  0.54279  2.28217 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              0.745755   0.359568   2.074 0.040129 *  
## Melanco_SMPIMelancolico -0.744325   0.208026  -3.578 0.000493 ***
## Intentossuicidio        -0.042617   0.029701  -1.435 0.153819    
## Psi_actualesSi           0.148143   0.288650   0.513 0.608700    
## MADRS                   -0.023358   0.011443  -2.041 0.043328 *  
## TBPTBP                  -0.216866   0.171993  -1.261 0.209692    
## IFD_B_EA2               -0.285788   0.191330  -1.494 0.137776    
## IFD_B_EA3               -0.614660   0.304731  -2.017 0.045831 *  
## IFD_B_EA4               -0.009464   0.378184  -0.025 0.980076    
## IFD_B_EA5               -1.300083   0.478203  -2.719 0.007486 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9151 on 125 degrees of freedom
## Multiple R-squared:  0.2662, Adjusted R-squared:  0.2134 
## F-statistic: 5.038 on 9 and 125 DF,  p-value: 8.479e-06

Anova(modelo_sin_TMTA, type = "III")

Observamos un aumento en el coeficiente beta estimado para la variable explicativa principal Melanco_SMPI que pasa de -0.49 a -0.74

Por otro lado, para comprender el funcionamiento de la variable TMTA sobre el resto de las variables de modelo final propuesto para FE, ajustamos un modelo con esas variables y TMTA como variable respuesta.

modelo_TMTA <- lm(TMTA ~ Melanco_SMPI + Intentossuicidio + Psi_actuales +
                     MADRS + TBP + IFD_B_EA,
                   data = melanco2)
summary(modelo_TMTA)

## 
## Call:
## lm(formula = TMTA ~ Melanco_SMPI + Intentossuicidio + Psi_actuales + 
##     MADRS + TBP + IFD_B_EA, data = melanco2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7274 -0.6680  0.0744  0.6422  2.7500 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)              1.024334   0.430297   2.381   0.0188 *
## Melanco_SMPIMelancolico -0.495726   0.248946  -1.991   0.0486 *
## Intentossuicidio         0.005453   0.035543   0.153   0.8783  
## Psi_actualesSi          -0.171343   0.345429  -0.496   0.6207  
## MADRS                   -0.022773   0.013694  -1.663   0.0988 .
## TBPTBP                   0.021220   0.205825   0.103   0.9180  
## IFD_B_EA2               -0.299412   0.228966  -1.308   0.1934  
## IFD_B_EA3                0.104902   0.364673   0.288   0.7741  
## IFD_B_EA4               -0.047263   0.452575  -0.104   0.9170  
## IFD_B_EA5               -1.089388   0.572268  -1.904   0.0593 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.095 on 125 degrees of freedom
## Multiple R-squared:  0.1387, Adjusted R-squared:  0.07674 
## F-statistic: 2.238 on 9 and 125 DF,  p-value: 0.02368

Anova(modelo_TMTA, type = "III")

Observamos que las variables explicativas para FE no resultan significativas para explicar TMTA, excepto para Melanco_SMPI (beta estimado = -0.49, p = 0.049). El modelo para TMTA obtiene un R² ajustado de 0.08

Calculamos la diferencia de medias en el desempeño en FE entre grupos de pacientes deprimidos, al incluir o no en el modelo a la variable TMT-A

comp_modelo_final<-emmeans (modelo_final, pairwise ~ Melanco_SMPI)
comp_modelo_final

## $emmeans
##  Melanco_SMPI   emmean    SE  df lower.CL upper.CL
##  No melancolico -0.498 0.156 124   -0.806   -0.189
##  Melancolico    -0.984 0.165 124   -1.310   -0.658
## 
## Results are averaged over the levels of: Psi_actuales, TBP, IFD_B_EA 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                     estimate    SE  df t.ratio p.value
##  No melancolico - Melancolico    0.487 0.166 124   2.930  0.0040
## 
## Results are averaged over the levels of: Psi_actuales, TBP, IFD_B_EA

comp_modelo_sin_TMTA<-emmeans(modelo_sin_TMTA, pairwise ~ Melanco_SMPI)
comp_modelo_sin_TMTA

## $emmeans
##  Melanco_SMPI   emmean    SE  df lower.CL upper.CL
##  No melancolico -0.546 0.198 125   -0.939   -0.154
##  Melancolico    -1.291 0.205 125   -1.696   -0.885
## 
## Results are averaged over the levels of: Psi_actuales, TBP, IFD_B_EA 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                     estimate    SE  df t.ratio p.value
##  No melancolico - Melancolico    0.744 0.208 125   3.578  0.0005
## 
## Results are averaged over the levels of: Psi_actuales, TBP, IFD_B_EA

Cuando el modelo no incluye a la variable TMT-A (medida de velocidad de procesamiento), la diferencia de medias en el desempeño en FE entre pacientes melancólicos y no melancólicos es mayor. Sin embargo, el cambio no es tan marcado, pasa de ser 0.49 a 0.74.

8.2 Visualización de resultados del modelo final para FE

Modelos	Variables	R2ajus
modelo1	Todas	.507
modelo_automatico	TBP +Melanco_SMPI + Intentossuicidio + Psi_actuales+ MADRS + IFD_B_BZD + IFD_B_AD + IFD_B_EA + IFD_B_AP + TMTA	.512
modelo_final	Melanco_SMPI + Intentossuicidio + Psi_actuales + MADRS + TBP + IFD_B_EA + TMTA	.514

8.2.1 Gráficos ggpredict

library(ggeffects)
fig<-ggpredict(modelo_final)
plot(fig, add.data = T)

## $Melanco_SMPI

## 
## $Intentossuicidio

## 
## $Psi_actuales

## 
## $MADRS

## 
## $TBP

## 
## $IFD_B_EA

## 
## $TMTA

8.2.2 Gráfico final

comparaciones<-emmeans(modelo_final, pairwise ~ Melanco_SMPI)
estimaciones<- as.data.frame(comparaciones$emmeans)

ggplot(estimaciones, aes(x=Melanco_SMPI, y=emmean), fill=Melanco_SMPI) + 
  labs(x="Tipo de depresión", y="Funciones ejecutivas") + 
  geom_errorbar(aes(ymin=emmean-SE, ymax=emmean+SE), color="darkblue", width=0.1)+
  geom_point(shape=15, size=4, color="darkblue")+  
  ggtitle("Rendimiento en funciones ejecutivas según tipo de depresión", "Media (EE)")+ 
  theme_light()

TP Taller: Melancolía y funciones ejecutivas. Versión 22/10

Grupo 5: Marina Valerio, Julieta Lomastro, Alejandra Vega y Samanta Leiva

28/10/2022