#Preparación de datos

Cargamos los datos iniciales y exploramos por encima las variables:

df=read_spss("Estudio Piloto Iria.sav")

Las mediciones repetidas conviene tenerlas en un formato alternativo, (largo en lugar de ancho):

dfLong = df %>% gather(Variable,Valor, -PACIENTE,-Brazo,-Edadcirugía,-Empleo, -TécnicaquirLaparosccir.abierta,-PSA) %>%
           separate(Variable,c("Prefijo","Tiempo")) %>%
           mutate(Tiempo=as.integer(Tiempo)) %>%
            spread(Prefijo,Valor) %>%
  mutate(Grupo=factor(Brazo,labels=c("Exp","Ctrl")),
         TiempoFactor=factor(Tiempo)) %>% 
  as_tibble()

 dfLong %>% head() %>% knitr::kable()
PACIENTE Brazo Edadcirugía Empleo TécnicaquirLaparosccir.abierta PSA Tiempo ICQSF PT1h PT24h Grupo TiempoFactor
0 2 65 2 1 4.1 0 13 13 281 Ctrl 0
0 2 65 2 1 4.1 1 11 2 20 Ctrl 1
0 2 65 2 1 4.1 2 13 2 17 Ctrl 2
0 2 65 2 1 4.1 3 13 3 14 Ctrl 3
0 2 65 2 1 4.1 6 11 3 9 Ctrl 6
5 2 63 3 1 4.4 0 13 10 224 Ctrl 0 
lasVariables=c("ICQSF","PT1h","PT24h")

Las variables con las que vamos a trabajar son entonces: ICQSF, PT1h, PT24h

Seguramente van a ser de interés los cambios con respecto a Baseline (Tiempo==0). Vamos a crear esa base de datos:

dfBaseline=dfLong %>% select(PACIENTE,Grupo,lasVariables) %>% group_by(PACIENTE,Grupo) %>% summarise_all(first)
names(dfBaseline)[-c(1,2)]=str_c(names(dfBaseline[-c(1,2)]),"___baseline")

dfDif=dfLong%>% select(PACIENTE,Grupo,Tiempo,lasVariables) %>% left_join(dfBaseline, by=c("PACIENTE","Grupo")) 

for (v in lasVariables) dfDif[[v]]=dfDif[[v]]-dfDif[[str_c(v,"___baseline")]]
dfDif=dfDif %>% filter(Tiempo >0) %>% select(-ends_with("___baseline")) %>%
  mutate(TiempoFactor=factor(Tiempo))

dfDif %>% head() %>% knitr::kable()
PACIENTE Grupo Tiempo ICQSF PT1h PT24h TiempoFactor
0 Ctrl 1 -2 -11 -261 1
0 Ctrl 2 0 -11 -264 2
0 Ctrl 3 0 -10 -267 3
0 Ctrl 6 -2 -10 -272 6
5 Ctrl 1 -1 3 -171 1
5 Ctrl 2 5 -3 -199 2

Vamos a trabajar de una en una con cada una de ellas

vTrabajo="ICQSF"

1 Análisis de ICQSF

Empezamos con una inspección visual:

ggplot(dfLong, aes(x = ICQSF)) + 
  geom_histogram()+ facet_wrap(~Tiempo)+
  theme_classic()

ggplot(dfLong, aes(x = Tiempo, y=ICQSF,group=Tiempo)) + 
  geom_boxplot()+
  theme_classic()

A continuación, los cambios:

ggplot(dfDif, aes(x = ICQSF)) + 
  geom_histogram()+ facet_wrap(~Tiempo)+
  theme_classic()

ggplot(dfDif, aes(x = Tiempo, y=ICQSF,group=Tiempo)) + 
  geom_boxplot()+
  theme_classic()

1.1 ICQSF comparado por Grupos

ggplot(dfLong, aes(x = ICQSF, y = as.factor(Tiempo),color=Grupo,fill=Grupo)) + 
  geom_density_ridges(aes(point_color=Grupo),alpha=0.25,scale=1.1, jittered_points=T, point_size=0.8 , point_alpha=0.5)+theme_classic()

1.2 Cambios en ICQSF comparado por Grupos

ggplot(dfDif, aes(x = ICQSF, y = as.factor(Tiempo),color=Grupo,fill=Grupo)) + 
  geom_density_ridges(aes(point_color=Grupo),alpha=0.25,scale=1.1, jittered_points=T, point_size=0.8,point_alpha=0.5)+theme_classic()

vTrabajo="PT1h"

2 Análisis de PT1h

Empezamos con una inspección visual:

  ggplot(dfLong, aes(x = PT1h)) + 
  geom_histogram()+ facet_wrap(~Tiempo)+
  theme_classic()

ggplot(dfLong, aes(x = Tiempo, y=PT1h,group=Tiempo)) + 
  geom_boxplot()+
  theme_classic()

A continuación, los cambios:

  ggplot(dfDif, aes(x = PT1h)) + 
  geom_histogram()+ facet_wrap(~Tiempo)+
  theme_classic()

ggplot(dfDif, aes(x = Tiempo, y=PT1h,group=Tiempo)) + 
  geom_boxplot()+
  theme_classic()

2.1 PT1h comparado por Grupos

ggplot(dfLong, aes(x = PT1h, y = as.factor(Tiempo),color=Grupo,fill=Grupo)) + 
  geom_density_ridges(aes(point_color=Grupo),alpha=0.25,scale=1.1, jittered_points=T, point_size=0.8 , point_alpha=0.5)+theme_classic()

2.2 Cambios en PT1h comparado por Grupos

ggplot(dfDif, aes(x = PT1h, y = as.factor(Tiempo),color=Grupo,fill=Grupo)) + 
  geom_density_ridges(aes(point_color=Grupo),alpha=0.25,scale=1.1, jittered_points=T, point_size=0.8,point_alpha=0.5)+theme_classic()

vTrabajo="PT24h"

3 Análisis de PT24h

Empezamos con una inspección visual:

ggplot(dfLong, aes(x = PT24h)) + 
  geom_histogram()+ facet_wrap(~Tiempo)+
  theme_classic()

ggplot(dfLong, aes(x = Tiempo, y=PT24h,group=Tiempo)) + 
  geom_boxplot()+
  theme_classic()

A continuación, los cambios:

ggplot(dfDif, aes(x = PT24h)) + 
  geom_histogram()+ facet_wrap(~Tiempo)+
  theme_classic()

ggplot(dfDif, aes(x = Tiempo, y=PT24h,group=Tiempo)) + 
  geom_boxplot()+
  theme_classic()

3.1 PT24h comparado por Grupos

ggplot(dfLong, aes(x = PT24h, y = as.factor(Tiempo),color=Grupo,fill=Grupo)) + 
  geom_density_ridges(aes(point_color=Grupo),alpha=0.25,scale=1.1, jittered_points=T, point_size=0.8 , point_alpha=0.5)+theme_classic()

3.2 Cambios en PT24h comparado por Grupos

ggplot(dfDif, aes(x = PT24h, y = as.factor(Tiempo),color=Grupo,fill=Grupo)) + 
  geom_density_ridges(aes(point_color=Grupo),alpha=0.25,scale=1.1, jittered_points=T, point_size=0.8,point_alpha=0.5)+theme_classic()

#Contrastando las diferencias por Grupo. En primer usar usamos lo valores crudos en cada tiempo:

generaTablatTestPorGrupo_Tiempo(dfLong,"Grupo","TiempoFactor", lasVariables, columnas = c("n", "mediaet", "p.t", "ci95", "p.w")) %>%
  knitr::kable(booktabs =T) %>%
    kable_styling(full_width = FALSE,latex_options = c("HOLD_position"),font_size = 9) %>%
  collapse_rows(columns=c(1), latex_hline = c("major"))
Variable Tiempo n.Exp mediaet.Exp n.Ctrl mediaet.Ctrl p.t ci95 p.w
ICQSF 0 25 13.48±0.80 22 15.36±0.70 0.075. 1.88[-0.20,3.97] 0.079.
1 25 12.32±0.66 22 13.77±0.68 0.124 1.45[-0.42,3.32] 0.102
2 25 9.28±0.86 21 12.48±0.93 0.014* 3.20[0.69,5.70] 0.011*
3 25 5.68±0.86 20 12.20±0.77 <0.001* 6.52[4.24,8.80] <0.001*
6 23 3.87±0.84 18 9.94±1.12 <0.001* 6.07[3.31,8.84] <0.001*
PT1h 0 25 72.48±19.24 22 61.09±15.64 0.641 -11.39[-60.29,37.51] 0.741
1 25 26.76±6.43 22 56.00±16.26 0.098. 29.24[-5.79,64.27] 0.162
2 25 13.12±4.10 21 51.67±15.77 0.024* 38.55[5.62,71.47] 0.069.
3 25 4.64±1.83 20 35.30±9.91 0.005* 30.66[10.19,51.13] <0.001*
6 23 0.70±0.35 18 19.50±7.34 0.017* 18.80[3.75,33.86] <0.001*
PT24h 0 25 465.48±99.23 22 443.91±93.16 0.872 -21.57[-289.81,246.67] 0.983
1 25 258.60±66.65 22 352.64±111.65 0.464 94.04[-164.11,352.18] 0.565
2 25 128.64±38.98 21 276.38±81.28 0.104 147.74[-32.30,327.78] 0.242
3 25 27.32±11.69 20 196.70±65.40 0.016* 169.38[34.38,304.38] 0.003*
6 23 4.00±1.50 18 107.78±43.00 0.024* 103.78[15.58,191.97] <0.001*

Ahora estudiamos los cambios con respecto a baseline

En primer usar usamos lo valores crudos en cada tiempo:

generaTablatTestPorGrupo_Tiempo(dfDif,"Grupo","TiempoFactor", lasVariables, columnas = c("n", "mediaet", "p.intra", "p.t", "ci95", "p.w")) %>%
  knitr::kable(booktabs =T) %>%
  kable_styling(full_width = FALSE,latex_options = c("HOLD_position"),font_size = 8) %>%
  collapse_rows(columns=c(1), latex_hline = c("major"))
Variable Tiempo n.Exp mediaet.Exp p.intra.Exp n.Ctrl mediaet.Ctrl p.intra.Ctrl p.t ci95 p.w
ICQSF 1 25 -1.16±0.73 0.124 22 -1.05±0.60 0.098. 0.902 0.11[-1.75,1.98] 0.889
2 25 -4.20±0.89 <0.001* 21 -2.38±0.97 0.023* 0.164 1.82[-0.77,4.41] 0.187
3 25 -7.80±1.00 <0.001* 20 -2.45±0.89 0.012* <0.001* 5.35[2.72,7.98] <0.001*
6 23 -9.39±1.03 <0.001* 18 -4.61±1.16 <0.001* 0.003* 4.78[1.71,7.85] 0.010*
PT1h 1 25 -45.72±17.99 0.018* 22 -18.09±12.98 0.178 0.210 27.63[-16.17,71.43] 0.159
2 25 -59.36±17.97 0.003* 21 -25.57±13.02 0.064. 0.127 33.79[-10.03,77.61] 0.186
3 25 -67.84±18.86 0.001* 20 -35.90±12.86 0.012* 0.160 31.94[-13.17,77.05] 0.458
6 23 -68.48±20.82 0.003* 18 -53.50±12.18 <0.001* 0.529 14.98[-32.88,62.84] 0.703
PT24h 1 25 -206.88±69.73 0.007* 22 -123.41±76.12 0.120 0.413 83.47[-120.05,287.00] 0.757
2 25 -336.84±76.32 <0.001* 21 -220.67±62.69 0.002* 0.236 116.17[-78.61,310.95] 0.724
3 25 -438.16±98.07 <0.001* 20 -262.30±77.97 0.003* 0.158 175.86[-71.30,423.02] 0.732
6 23 -412.17±93.78 <0.001* 18 -317.11±57.22 <0.001* 0.381 95.06[-122.68,312.80] 0.937

4 Calculos de potencia

library(pwr)



pwr.t2n.test(n1=23,n2=18,d=4*18.80/(3.75-33.86))
## 
##      t test power calculation 
## 
##              n1 = 23
##              n2 = 18
##               d = 2.5
##       sig.level = 0.05
##           power = 1
##     alternative = two.sided
resumen=dfLong %>% filter(Tiempo==6) %>% select(Grupo,ICQSF,PT1h,PT24h) %>% pivot_longer(cols = c(ICQSF,PT1h,PT24h)) %>% group_by(Grupo,name) %>% summarise(mean=mean(value,na.rm=T),sd=sd(value,na.rm=T),n=sum(!is.na(value)))


cohen=resumen %>% group_by(name) %>% summarise(Sp=sqrt(first(sd)^2+last(sd)^2),dif=abs(first(mean)-last(mean))) %>% mutate(d=dif/Sp)

Cálculos de potencia post-hoc

pwr.t2n.test(n1=23,n2=18,d=cohen %>% filter(name=="ICQSF") %>% .[["d"]])
## 
##      t test power calculation 
## 
##              n1 = 23
##              n2 = 18
##               d = 1
##       sig.level = 0.05
##           power = 0.87
##     alternative = two.sided
pwr.t2n.test(n1=23,n2=18,d=cohen %>% filter(name=="PT1h") %>% .[["d"]])
## 
##      t test power calculation 
## 
##              n1 = 23
##              n2 = 18
##               d = 0.62
##       sig.level = 0.05
##           power = 0.49
##     alternative = two.sided
pwr.t2n.test(n1=23,n2=18,d=cohen %>% filter(name=="PT24h") %>% .[["d"]])
## 
##      t test power calculation 
## 
##              n1 = 23
##              n2 = 18
##               d = 0.58
##       sig.level = 0.05
##           power = 0.44
##     alternative = two.sided

La potencia post-hoc para ICQSF a los 6 meses ha sido del 87%,que se suele considerar como una muestra de buen tamaño.

Para PT24h y PT1h la cosa no resulta tan buena (44% y 49%). En general para resultados como los obtenidos, en estudios futuros debería considerarse muestras más grandes.