knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
library(data.table)
## Warning: package 'data.table' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'data.table'
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
library(haven)

library(geepack)
## Warning: package 'geepack' was built under R version 4.4.3
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.4.3
## corrplot 0.95 loaded
library(psych)
## 
## Adjuntando el paquete: 'psych'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
#library(sjPlot) opcional para ver resultados = tab_model(modelo)
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
options(scipen=999, digits=3)
CD4 <- read_dta("C:/Users/Administrador/Downloads/CD4.dta")
#View(CD4_Ejercicio1)

colSums(is.na(CD4)) ### Valores Faltantes 
##        id      time       cd4       Age    gender       sex Treatment 
##         0         0        42         0         0         0         0

La base esta en formato “LONG” Los datos se muestran en formato “long”, donde cada fila es una medición y time es el momento de la medición. El id se repite para cada individuo de acuerdo al n de mediciones que tuvo. 1)Evaluar el cambio de CD4 en el tiempo ajustando por edad y sexo como posibles confundidores 2) Evaluar el cambio de CD4 en función del tiempo EN HOMBRES Y MUJERES 3) Evaluar el cambio de CD4 en función del tiempo en los grupos de TRATAMIENTO

CD4$Treatment <- as.numeric(CD4$Treatment)  
# LINK: https://tidyr.tidyverse.org/articles/pivot.html#wider
# Para ver las correlaciones, debemos pasar lods datos a formato "WIDE"

CD4_WIDE <- CD4 %>% 
  pivot_wider(names_from = time, ## selecciono la variable repetida
              names_prefix = "Tiempo", ## agrego el prefijo para el nuevo nombre
              values_from = cd4)  ### selecciono de donde va a tomar los valores

# AL PASAR DE LONG A WIDE se  ELIMINA LA VARIABLE TIME Y SE CREAN LAS 3 
# MEDIDAS REPETIDAS DE CD4 CON LA INFORMACION DE TIME, Tiempo0, Tiempo6 Y Tiempo12

CD4_WIDE[1:10,-3] 
## # A tibble: 10 × 7
##    id             Age   sex Treatment Tiempo12 Tiempo6 Tiempo0
##    <dbl+lbl>    <dbl> <dbl>     <dbl>    <dbl>   <dbl>   <dbl>
##  1 18 [DIC-107]    23     1         0      529     281     364
##  2 47 [JB-11]      30     1         0      506     605     311
##  3 11 [BKC-32]     33     1         0      735     780     554
##  4 71 [MLL-61*]    47     1         0      670     638     468
##  5 69 [ML-37]      44     1         0      365     487     281
##  6 65 [MDL-81]     38     0         0     1024     850     704
##  7 14 [CCG-45]     35     1         0      613     524     377
##  8 75 [NEC-119]    31     1         0      729     590     521
##  9 20 [DK-04]      34     1         0      682     777     652
## 10 62 [LFD-39]     26     1         0      434     372     607
View(CD4_WIDE)
head(CD4[order(CD4$id, decreasing = F),]) %>% kable() 
id time cd4 Age gender sex Treatment
1 0 311 62 M 1 1
1 6 473 62 M 1 1
1 12 408 62 M 1 1
2 0 436 25 M 1 0
2 6 442 25 M 1 0
2 12 643 25 M 1 0
view(CD4)
corrplot.mixed(cor(CD4_WIDE[,6:8], # selecciono las variables dependientes correlacionadas
                   use = "pairwise.complete"), # correlacion
               lower = "number", # panel inferior = número
               upper = "circle", # panel superior = Círculos
         title = "Correlacion CD4") # título

CD4_Long <- CD4_WIDE %>%  # selecciono base wide
  pivot_longer(
    cols = starts_with("Tiempo"), #  columnas que empiezan con tiempo (0, 6 ,12)
    names_to = "tiempo_nuevo", #  Nueba variable que va a contener tiempo(0, 6, 12)
    names_transform = as.integer, # transformo tiempo a integer-numerica
    names_prefix = "Tiempo",  # saco la palabra tiempo para que no quede tiempo6 como valor
    values_to = c("cd4_nuevo") # los valores de cd4 de cada tiempo a la nueva variable 
  )

head(CD4_Long)
## # A tibble: 6 × 7
##   id             Age gender   sex Treatment tiempo_nuevo cd4_nuevo
##   <dbl+lbl>    <dbl> <chr>  <dbl>     <dbl>        <int>     <dbl>
## 1 18 [DIC-107]    23 M          1         0           12       529
## 2 18 [DIC-107]    23 M          1         0            6       281
## 3 18 [DIC-107]    23 M          1         0            0       364
## 4 47 [JB-11]      30 M          1         0           12       506
## 5 47 [JB-11]      30 M          1         0            6       605
## 6 47 [JB-11]      30 M          1         0            0       311
#modelo de regresión lineal si las variables fueran independientes
mod_lineal <- lm(cd4~time +Age + sex, data = CD4)  

summary(mod_lineal) 
## 
## Call:
## lm(formula = cd4 ~ time + Age + sex, data = CD4)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -531.2 -149.0  -18.7  123.7  748.3 
## 
## Coefficients:
##             Estimate Std. Error t value      Pr(>|t|)    
## (Intercept)  479.846     71.738    6.69 0.00000000013 ***
## time          17.224      2.703    6.37 0.00000000083 ***
## Age            0.857      1.529    0.56          0.58    
## sex          -13.929     55.412   -0.25          0.80    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 217 on 263 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.135,  Adjusted R-squared:  0.125 
## F-statistic: 13.7 on 3 and 263 DF,  p-value: 0.000000026
#modelo asumiendo dependencia
mod_gee <- geeglm(cd4 ~ time + Age + sex ,   # y~ covariables
                     data = CD4,   # datos
                     id = CD4$id,   # id= identificador del grupo al que corresponden las                                                  medidas repetidas
                    family = gaussian,         # modelo lineal (existen otras: poisson,gaussian, etc)
                    corstr = "unstructured") # tipo de correlación ("exchangeable", "independence" ,                                                "ar1", etc)


tidy(mod_gee, conf.int = T)
## # A tibble: 4 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)  481.        63.6    57.3    3.77e-14   357.      606.  
## 2 time          17.6        2.81   39.2    3.78e-10    12.1      23.1 
## 3 Age            0.734      1.43    0.265  6.07e- 1    -2.06      3.53
## 4 sex          -13.1       48.6     0.0722 7.88e- 1  -108.       82.3

¿COMO SE INTERPRETAN LOS COEFICIENTES? ¿QUE PASA CON LOS BETA Y LOS ERRORES STANDARD EN UNO Y OTRO MODELO? VER EFECTO DEL TIEMPO: ¿CUAL ES LA MAGNITUD DE CAMBIO DE CD4 POR MES? Y EL DE GRUPOS ¿EXISTE DIFERENCIA BASAL DE CD4 POR SEXO?

El número de cd4 en tiempo cero cuando el sexo es femenino es 581.4 El número de cd4 aumenta en 17.6 por cada mes El número de cd4 no se modifica en forma estadísticamente significativa según la edad El número de cd4 basal no se modifica en forma estadísticamente significativa según el sexo Los errores disminuyen en el modelo GEE

#GRÁFICO DE PREDICCIÓN LINEAL DEL EFECTO DEL TIEMPO SOBRE CAMBIO DE CD4 CON IC95%

#Observados

ggplot(CD4, aes(time,cd4 )) + geom_point()+
 geom_smooth(method = "lm", se = TRUE) + 
  theme_classic()+ 
 labs(x = "Tiempo", y = "CD4", title = "Efectos marginales") 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 42 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 42 rows containing missing values or values outside the scale range
## (`geom_point()`).

#GRÁFICO DE PREDICCIÓN LINEAL DEL EFECTO DEL TIEMPO SOBRE CAMBIO DE CD4 

# Comparo predichos

datospredichosGEE <- data.frame(na.omit(CD4), ## sacamos datos faltantes 
                                                         # o no podemos graficar
                                phat = fitted(mod_gee))

ggplot(datospredichosGEE, aes(time, phat)) + 
 geom_smooth(method = "lm") +
ylim (20, 1400) +
  theme_classic()+ 
 labs(x = "Tiempo", y = "CD4", title = "Efectos marginales") 
## `geom_smooth()` using formula = 'y ~ x'

Evaluacion entre sexo y tiempo

mod_gee2 <- geeglm(cd4 ~ time*sex  ,
                     data = CD4,
                     id = CD4$id,
                     corstr = "unstructured")
tidy(mod_gee2, conf.int = T) 
## # A tibble: 4 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)   486.       75.8    41.2    1.40e-10   338.       635. 
## 2 time           20.7       9.95    4.32   3.78e- 2     1.17      40.2
## 3 sex            10.1      78.5     0.0166 8.98e- 1  -144.       164. 
## 4 time:sex       -3.24     10.4     0.0973 7.55e- 1   -23.6       17.1

¿EL CAMBIO DE CD4 EN EL TIEMPO ES DIFERENTE EN HOMBRES Y MUJERES? no

#Evaluación del efecto del tratamiento

#Analisis exploratorio: media de CD4 en cada medición por grupo de tratamiento

CD4_WIDE %>% group_by(Treatment) %>% summarise(T0 = mean(Tiempo0, na.rm = T), 
                                               T6 = mean(Tiempo6, na.rm = T), 
                                               T12 = mean(Tiempo12, na.rm = T))
## # A tibble: 2 × 4
##   Treatment    T0    T6   T12
##       <dbl> <dbl> <dbl> <dbl>
## 1         0  535.  601.  648.
## 2         1  468.  631.  730.
mod_gee2 <- geeglm(cd4 ~ time +Treatment+ sex  ,
                     data = CD4,
                     id = CD4$id,
                     corstr = "unstructured")
tidy(mod_gee2, conf.int = T) 
## # A tibble: 4 × 7
##   term        estimate std.error statistic  p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 (Intercept)   495.       54.1    83.7    0           389.      601. 
## 2 time           17.7       2.84   39.0    4.15e-10     12.2      23.3
## 3 Treatment      11.9      26.2     0.206  6.50e- 1    -39.5      63.3
## 4 sex            -7.77     48.7     0.0254 8.73e- 1   -103.       87.7

¿MIRANDO LOS VALORES A LOS 0 6 Y 12 MESES, LE PARECE QUE EXISTEN DIFERENCIAS ENTRE LOS GRUPOS? Parecería que sí

MODELO MARGINAL (PROMEDIO POBLACIONAL) DE REGRESION LINEAL ASUMIENDO DEPENDENCIA

VER EFECTO DEL TIEMPO (MAGNITUD DEL CAMBIO POR MES) Y EL DEL TRATAMIENTO (¿EXISTE DIFERENCIA BASAL DE CD4 ENTRE LOS GRUPOS?)No existe diferencia basal en los grupos de tratamiento

#gráfico de trayectoria de tratados y no tratados
mod_gee4 <- geeglm(cd4 ~ time*Treatment + sex ,
                     data = CD4, 
                     id = CD4$id,
                     corstr = "unstructured")
tidy(mod_gee4, conf.int = T) 
## # A tibble: 5 × 7
##   term           estimate std.error statistic p.value conf.low conf.high
##   <chr>             <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 (Intercept)      544.       53.5   104.     0         440.       649. 
## 2 time               9.42      3.65    6.64   0.00998     2.25      16.6
## 3 Treatment        -63.6      37.9     2.81   0.0935   -138.        10.7
## 4 sex               -7.06     49.1     0.0207 0.886    -103.        89.1
## 5 time:Treatment    13.2       5.29    6.23   0.0126      2.83      23.6
#GRAFICO DEL CAMBIO DE CD4 en tratados y no tratados

predichosTTO <- data.frame(na.omit(CD4), TTO = fitted(mod_gee4))

ggplot(predichosTTO, aes(time, TTO, colour = factor(Treatment))) + 
 geom_smooth(method = "lm") + theme_classic()+ ggtitle("Efectos marginales") +
 labs(x = "Tiempo", y = "CD4") 
## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

¿EL CAMBIO DE CD4 EN EL TIEMPO ES DIFERENTE EN LOS GRUPOS DE TRATAMIENTO? ¿CUÁL ES LA MAGNITUD DE CAMBIO POR MES EN CADA GRUPO? ¿CUÁL ES EL VALOR ESPERADO DE CD4 AL AÑO EN CADA GRUPO?

Es diferente. En el grupo sin tratamiento, los cd4 aumentan 9.42 cd4 por mes. en el grupo con tratamiento, aumentan 13.19+9.42=22.61