Aproximación aplicada y autodidacta al análisis de experimentos: ANOVA mediante el programa R

---

Cita: Ramos-Álvarez, M.M. (2017).Aproximación aplicada y autodidacta al análisis de experimentos: ANOVA mediante el programa R. Recuperado 14 diciembre 2017, desde http://www4.ujaen.es/~mramos/ADMmRa/ANOVAEnRMmRa.pdf

---

Ejemplo extraído de: Maxwell, S.E. & Delaney, H.D. (2004, 2Ed). Designing Experiments and Analyzing Data.
Although different mood states have, of course, always been of interest to clinicians, recent years have seen a profusion of studies attempting to manipulate mood states in controlled laboratory studies. In such induced-mood research, participants typically are randomly assigned to one of three groups: a depressed-mood induction, a neutral-mood induction, or an elated-mood induction. One study (Pruitt, 1988) used selected videoclips from several movies and public television programs as the mood-induction treatments. After viewing the video for her assigned condition, each participant was asked to indicate her mood on various scales. In addition, each subject was herself videotaped, and her facial expressions of emotion were rated on a scale of 1 to 7 (1 indicating sad; 4, neutral; and 7, happy) by an assistant who viewed the videotapes but was kept “blind” regarding the subjects’ assigned conditions. Table 3.3 shows representative data7 of these Global Affect Ratings for 10 observations per group, along with the means and standard deviations for the groups. These are the data displayed in Figure 3.1 on page 68.
Se trata de un Diseño Unifactorial EntreGrupos (DUEG) puesto que los 3 grupos corresponden a muestras independientes.
Fichero de datos: http://www4.ujaen.es/~mramos/ADMmRa/MaxwelDelaney04Table3_3.csv

Código del tutorial en R: http://www4.ujaen.es/~mramos/ADMmRa/CodeANOVAEnR.R

1 R como complemento de otros programas (e.g. SPSS)

Si únicamente vamos a utilizar R como complemento de otros programas, sería algo así como una calculadora muy potente a la que vamos añadiendo funcionalidad al indicarle las librerías específicas que nos servirán para llevar a cabo aquellos análisis específicos que nos interesan. La fuente de referencia básica para encontrar las librerías será la página oficial del programa R (https://cloud.r-project.org/web/packages/).

# Por ejemplo, para llevar a cabo estimaciones del Tamaño del efecto
# y de la Potencia cargamos las librerías de la siguiente manera:
require("compute.es")
require("pwr")
# A partir de aquí incluiremos los valores concretos de los estadísticos
# que vayamos necesitando, que obtendremos de las pantallas de salida del
# programa SPSS por ejemplo.
LosTam=c(10,10,10)
LasVarAOV<-c(23.333333,  0.962963)
# Hemos introducido los tamaños muestrales de los 3 grupos del Diseño,
# así como las Varianzas (ó Medias de Cuadrados) que se obtienen de la
# Tabla ANOVA correspondiente.
# La función "power.anova.test" de la librería {pwr} nos permite estimar
# la potencia estadística asociada a los datos del diseño.
power.anova.test(groups = 3, n = LosTam[1], between.var = LasVarAOV[1], 
                 within.var = LasVarAOV[2])

## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 3
##               n = 10
##     between.var = 23.33333
##      within.var = 0.962963
##       sig.level = 0.05
##           power = 1
## 
## NOTE: n is number in each group

Mejor intentaremos sacar todo el partido posible al Programa R

---

2 Preparativos

Primero cargamos las librerías,

# La que es específica del curso
CarpAct<-"http://www4.ujaen.es/~mramos/ADMmRa/"
FileAct<-"FuncApoyoADMmRa.R"
source(paste0(CarpAct,FileAct))
# Y todas aquellas que usaremos, obtenidas en su mayoría del
#servidor oficial de R Project:
initPkg(c("MBESS", "metafor", "plyr","data.table","car","effsize",
          "compute.es","bootES","pwr","ez","coin","psych","lsmeans",
          "lsr","Rmisc","gplots","ggplot2","sm","Hmisc","forestplot"))

Y el fichero de datos para trabajar:

CarpAct<-"http://www4.ujaen.es/~mramos/ADMmRa/"
# En la Variable CarpAct se especifica la ruta.
# En Windows (probablemente cambie el usuario aulas):
  # Carpeta descargas: CarpAct<-"C:/Users/aulas/Downloads/"
  # Carpeta Documentos: CarpAct<-"C:/Users/aulas/Documents/"
  # Escritotio: CarpAct<-"C:/Users/aulas/Desktop/"
FileAct<-"MaxwelDelaney04Table3_3.csv"
DataRaw<-read.table(paste0(CarpAct,FileAct),header=T,sep=";",dec=",")
Datos<-data.table(DataRaw) #Requiere la librería {data.table}

Nuestros datos están contenidos en la variable que hemos denominado Datos,
definida según el formato data.table, derivada a su vez del formato más importante de R, data.frame.
Para visualizar el contenido de la misma, escribimos el nombre de la variable y lo ejecutamos (Comando Run o Excecute),

Datos

##     IdSubj      Group DV
##  1:      1   Pleasant  6
##  2:      2   Pleasant  5
##  3:      3   Pleasant  4
##  4:      4   Pleasant  7
##  5:      5   Pleasant  7
##  6:      6   Pleasant  5
##  7:      7   Pleasant  5
##  8:      8   Pleasant  7
##  9:      9   Pleasant  7
## 10:     10   Pleasant  7
## 11:     11    Neutral  5
## 12:     12    Neutral  4
## 13:     13    Neutral  4
## 14:     14    Neutral  3
## 15:     15    Neutral  4
## 16:     16    Neutral  3
## 17:     17    Neutral  4
## 18:     18    Neutral  4
## 19:     19    Neutral  4
## 20:     20    Neutral  5
## 21:     21 Unpleasant  3
## 22:     22 Unpleasant  3
## 23:     23 Unpleasant  4
## 24:     24 Unpleasant  4
## 25:     25 Unpleasant  4
## 26:     26 Unpleasant  3
## 27:     27 Unpleasant  1
## 28:     28 Unpleasant  2
## 29:     29 Unpleasant  2
## 30:     30 Unpleasant  4
##     IdSubj      Group DV

---

3 Análisis convencionales del ANOVA

3.1 Base del ANOVA: Definir el Modelo

Modelo=as.formula(DV~ Group)
M.AOV=aov(Modelo,Datos)

---

3.2 Realizar el ANOVA

ResAOV<-summary(M.AOV)
ResAOV2<-anova(M.AOV)
ResLM<-summary.lm(M.AOV) # ANOVA como Modelo Lineal
# Visualizamos las Variables creadas para comprobar las diferencias
ResAOV

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Group        2  46.67  23.333   24.23 9.42e-07 ***
## Residuals   27  26.00   0.963                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ResAOV2

## Analysis of Variance Table
## 
## Response: DV
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Group      2 46.667  23.333  24.231 9.421e-07 ***
## Residuals 27 26.000   0.963                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ResLM

## 
## Call:
## aov(formula = Modelo, data = Datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##     -2     -1      0      1      1 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       4.0000     0.3103  12.890 4.75e-13 ***
## GroupPleasant     2.0000     0.4389   4.557   0.0001 ***
## GroupUnpleasant  -1.0000     0.4389  -2.279   0.0308 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9813 on 27 degrees of freedom
## Multiple R-squared:  0.6422, Adjusted R-squared:  0.6157 
## F-statistic: 24.23 on 2 and 27 DF,  p-value: 9.421e-07

# Extraemos los términos más importantes del Diseño y de la Tabla ANOVA
# Nota: Es conveniente familiarizarse con los contenidos, pulsando 
# sobre el nombre de cualquiera de las variables y ejecutándolo
LosNivBetw=levels(Datos[,Group])
NNivBetw=length(LosNivBetw) #Nº grupos
nBetw=Datos[, .N, by = Group] #Tamaños muestrales para cada grupo
LasMed<-Datos[,mean(DV),by=Group][[2]]
LosTam=Datos[,.N,by=Group][[2]]
LasVarAOV<-ResAOV2[, "Mean Sq", drop = FALSE][[1]]
VarErr=last(LasVarAOV)
LaF<-ResAOV2[, "F value", drop = FALSE][[1]][1]
laP<-ResAOV2[, "Pr(>F)", drop = FALSE][[1]][1]
LosGradosLib<-ResAOV2[, "Df", drop = FALSE][[1]]
NTT=Datos[,.N]

---

3.3 Análisis Post Hoc

# Una vez realizado el Análisis Global (o perspectiva Omnibus),
# normalmente las revistas científicas esperan encontrar el análisis Post Hoc ...
# Y en este diseño básico tipo "Entre", la prueba más popular es la de Tukey:
posthoc <- TukeyHSD(x=M.AOV)
posthoc

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Modelo, data = Datos)
## 
## $Group
##                     diff        lwr         upr     p adj
## Pleasant-Neutral       2  0.9118983  3.08810169 0.0002865
## Unpleasant-Neutral    -1 -2.0881017  0.08810169 0.0761816
## Unpleasant-Pleasant   -3 -4.0881017 -1.91189831 0.0000007

---

3.4 El análisis descriptivo mediante una representación gráfica

# No se nos había olvidado. En realidad, el análisis descriptivo es previo,
# pero es conveniente entender el modelo antes de nada y por esto 
# hemos alterado el orden de presentación.
# En R la representación de los promedios no corresponde por defecto
# a un simple diagrama de Barras, motivo por el cual tendremos que trabajar
# un poco en el programa para realizar lo que muchas revistas nos exigirán
# en este punto, aunque esté desfasado (Difícil, Seleccionar, Una, Unica & Cita; 2017) ...
# Previamente obtenemos los estadísticos descriptivos:
LosEstad <- Datos[,list(M=mean(DV),Sd=sd(DV)),by=Group]
# Para llevarlos entonces al gráfico correspondiente,
# Si la VarInd es categórica en origen:
Graph<-LosEstad[,barplot(M, names.arg= Group,ylim=c(0,10))];

# Si, por el contraio, la VarInd es cuantitativa en origen (que no es el caso
# de nuestro diseño pero lo haremos con fines didácticos, para no cambiar de
# supuesto de prácticas):
plotmeans(Datos$DV~Datos$Group) # Requiere la librería {gplots}

#Terminamos con la exploración de los estadísticos descriptivos
LosEstad

##         Group M        Sd
## 1:   Pleasant 6 1.1547005
## 2:    Neutral 4 0.6666667
## 3: Unpleasant 3 1.0540926

Sería recomendable que anotase la interpretación de los resultados del Análisis realizado

## Interpretación 1:  
##  
##  
##  
##  
##  
##  
##  
##  
##  
##  
###################################

---

4 Ampliación de los Análisis del ANOVA

Ahora vamos más allá de lo convencional …
El objetivo es comprender los datos de mi experimento en profundidad, bucear (diving) en los mismos.
Recuerda: Primero describir (entender) y entonces inferir (concluir).

---

4.1 Descripción de los datos

# Obtenemos Estadísticos Descriptivos robustos (e.g. Mediana) en comparación
# a los clásicos, más el Error Típico de la Media (ETM) como medida de
# variabilidad Error,
StatsMMv2(DatDT=Datos,"Group","DV") # Requiere la librería general {FuncMmRA.R}

##         Group DV.N DV.Mean DV.Md DV.ETM DV.CI_95
## 1:   Pleasant   10       6   6.5  0.365    0.826
## 2:    Neutral   10       4   4.0  0.211    0.477
## 3: Unpleasant   10       3   3.0  0.333    0.754

---

4.2 Análisis Exploratorio de los datos

El objetivo es llevar la estadística descriptiva a un enálisis exploratorio basado en la aproximación estadística Exploratory Data Analysis (EDA),

# Primero con un Diagrama de Cajas y Barbas
boxplot(Modelo,Datos)

# Alternativamente, si optamos por los gráficos más convencionales,
LosEstadMas <- Datos[,list(M=mean(DV),Sd=sd(DV),ETM=ETM(DV),CI=CI_MM(DV)),by=Group]
# Añadimos la barra de Error que corresponde al ETM, requiere la librería {Hmisc}
Graph<-LosEstadMas[,barplot(M, names.arg= Group,ylim=c(0,10))];
LosEstadMas[errbar(Graph,M,M+ETM,M-ETM, add=TRUE)];

## Null data.table (0 rows and 0 cols)

# Añadimos la barra de Error que corresponde al Intervalo Confidencial (CI), que es
# tal vez menos recomendable pues se trata de describir, pero es lo que nos pedirán
# las revistas de manera mayoritaria:
Graph<-LosEstadMas[,barplot(M, names.arg= Group,ylim=c(0,10))];
LosEstadMas[errbar(Graph,M,M+CI,M-CI, add=TRUE)];

## Null data.table (0 rows and 0 cols)

#De manera alternativa, pueden obtener las estimaciones CI
# a partir de la librería {Rmisc}:
summarySE(data=Datos,measurevar="DV",groupvars="Group")

##        Group  N DV        sd        se        ci
## 1    Neutral 10  4 0.6666667 0.2108185 0.4769046
## 2   Pleasant 10  6 1.1547005 0.3651484 0.8260230
## 3 Unpleasant 10  3 1.0540926 0.3333333 0.7540524

# Para explorar los supuestos nos valdremos de lo que se conoce como
# gráficos para el diagnóstico del Modelo:
# (1) un gráfico de densidad de Kernel,
Datos[,sm.density.compare(DV, Group)] #Requiere la librería  {sm}

## NULL

IndVar=Datos$Group
colfill<-c(2:(2+length(levels(IndVar))));legend(0,0.5, levels(IndVar),
        fill=colfill) #Le añade una leyenda a la izquierda

# Y además un (2) Gráfico Q-Q Normal para explorar el supuesto de Normalidad.
  # Opción 1: Todos los grupos juntos
  qplot(sample = DV, data = Datos, color=Group) #Requiere la librería  {ggplot2}

  #Opción 2: Un gráfico por cada Grupo separado (valiéndonos de un bucle)
    lapply(LosNivBetw, function(x) {
    PorGrupos <- Datos[Group==x,DV] # Para ir introduciendo los niveles
    q<-qqPlot(PorGrupos,dist= "norm",main =x) # Requiere la librería  {car}
  })

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL

Un artículo concreto nos servirá para ilustrar estas posibilidades:
Martín-Guerrero, T.L., Ramos-Álvarez, M.M., and Rosas, J.M. (2016). Payoff affects tasters’ decisions but it does not affect their sensitivity to basic tastes. Food Quality and Preference, 48, 11-16. Doi: http://dx.doi.org/10.1016/j.foodqual.2015.07.016

… El Mundo gráfico de R es impresionante, probar por ejemplo la librería {ggplot2}
Este es un buen momento para realizar de nuevo la interpretación de los resultados del Análisis realizado

## Interpretación 2:  
##  
##  
##  
##  
##  
##  
##  
##  
##  
##  
###################################

---

4.3 Estimaciones relacionadas con el estadístico de Reducción Proporcional del Error (RPE) y la Potencia

Antes de empezar, una visión interpretativa general de este concepto:

Estimamos ahora el Tamaño del Efecto y sus Intervalos Confidenciales asociados.

# Para obtener la variante más usual del ANOVA
etaSquared(M.AOV) #Requiere la librería {lsr}

##          eta.sq eta.sq.part
## Group 0.6422018   0.6422018

# Alternativa propia con más Estimadores, que requiere tener previamente 
# cargada la librería {FuncMmRA.R}
EfSizeAOV(M.AOV) 

##          eta.sq eta.sq.part Omega.part  R.Sq.adj
## Group 0.6422018   0.6422018  0.6076459 0.6156983

# Para añadir los IC del ANOVA global: ¡¡¡ Ojo, la probabilidad tiene que
# estar fijada en .90 !!! Se basa en las funciones de la librería {MBESS}
# que está en revisión pues arroja errores.
ci.pvaf(F.value=LaF, df.1=LosGradosLib[1], df.2=LosGradosLib[2], N=NTT,conf.level=.90)

## $Lower.Limit.Proportion.of.Variance.Accounted.for
## [1] 0.4090366
## 
## $Probability.Less.Lower.Limit
## [1] 0.05
## 
## $Upper.Limit.Proportion.of.Variance.Accounted.for
## [1] 0.734783
## 
## $Probability.Greater.Upper.Limit
## [1] 0.05
## 
## $Actual.Coverage
## [1] 0.9

# Finalmente, obtenemos la potencia estadística para nuestro estudio.
power.anova.test(groups = 3, n = LosTam[1], between.var = LasVarAOV[1],
                 within.var = LasVarAOV[2]) # Requiere la librería {pwr}

## 
##      Balanced one-way analysis of variance power calculation 
## 
##          groups = 3
##               n = 10
##     between.var = 23.33333
##      within.var = 0.962963
##       sig.level = 0.05
##           power = 1
## 
## NOTE: n is number in each group

La estimación de los Intervalos Confidenciales asociados a las Medidas del Tamaño del Efecto es lo más, se puede constatar incluyendo en el buscador de Google los términos “effect size confidence interval”

---

4.4 Análisis detallado de Contrastes

El Análisis Global de los datos proporciona una visión muy burda en los diseños con más de dos grupos. Usualmente nos haremos preguntas mucho más concretas, de manera que para dedicamos ahora unos minutos a e este concepto:

En R, nos apoyaremos en una librería muy completa, {lsmeans}, que nos servirá además para otros
diseños y otras variantes del ANOVA más complejos

---

4.4.1 Definir el Modelo de Contrastes

ModContr <- lsmeans(M.AOV, "Group")
plot(ModContr,xlab="DV") # Representación de los Intervalos Confidenciales

---

4.4.2 Análisis con enfoque planeado (o también llamado “a priori”)

MisContr=list(c1=c(1,0,-1), c2=c(-1/2,1,-1/2))
ResCon1<-contrast(ModContr,MisContr , by = NULL)
test(ResCon1)

##  contrast estimate        SE df t.ratio p.value
##  c1            1.0 0.4388537 27   2.279  0.0308
##  c2            2.5 0.3800585 27   6.578  <.0001

confint(ResCon1)

##  contrast estimate        SE df   lower.CL upper.CL
##  c1            1.0 0.4388537 27 0.09954653 1.900453
##  c2            2.5 0.3800585 27 1.72018442 3.279816
## 
## Confidence level used: 0.95

plot(ResCon1) 

# Recuerda que puedes añadir xlab="DV", para cambiar los nombres en los ejes

Ahora, un poco de teoría vendrá bien para entender las entrañas de los contrastes

---

4.4.3 Análisis de Tendencias

Si el diseño incluyese una variable independiente cuantitativa en origen,
entonces el análisis detallado se podría enfocar como si de un Modelo Lineal (y sus variantes)
se tratase.
En el ejemplo no sería adecuado este tipo de análisis puesto que la variable es categórica, no obstante aplicaremos la técnica para no tener que cargar un nuevo fichero de datos. Imaginemos por un momento que se trata de 3 tipos de imagen caracterizadas con niveles de agrado bajo, medio y elevado, según una escala de valencia de Lang)

ResCon2<-contrast(ModContr, "poly", name = "order")
test(ResCon2)

##  order     estimate        SE df t.ratio p.value
##  linear          -1 0.4388537 27  -2.279  0.0308
##  quadratic       -5 0.7601170 27  -6.578  <.0001

confint(ResCon2)

##  order     estimate        SE df  lower.CL    upper.CL
##  linear          -1 0.4388537 27 -1.900453 -0.09954653
##  quadratic       -5 0.7601170 27 -6.559631 -3.44036885
## 
## Confidence level used: 0.95

plot(ResCon2)

¿En qué consiste el Análisis de Tandencias?

---

4.4.4 Análisis A Posteriori (ó no planeado), poupularizado como tipo Post Hoc

Ahora el objetivo es desmenuzar con lupa todo el diseño, rastreando las posibles diferencias que surgen al comparar todos los niveles unos con otros …

ResCon3<-contrast(ModContr, "pairwise") #Por defecto asume Tukey
test(ResCon3)

##  contrast              estimate        SE df t.ratio p.value
##  Neutral - Pleasant          -2 0.4388537 27  -4.557  0.0003
##  Neutral - Unpleasant         1 0.4388537 27   2.279  0.0762
##  Pleasant - Unpleasant        3 0.4388537 27   6.836  <.0001
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

 # Para aplicar otro tipo de ajustes:
Adjust2<-test(ResCon3, side = "=",  adjust = "holm")
Adjust3<-test(ResCon3, side = "=",  adjust = "bonferroni")
Adjust4<-test(ResCon3, side = "=",  adjust = "fdr") # Muy recomnedable: Benjamini & Hochberg
confint(ResCon3)

##  contrast              estimate        SE df    lower.CL   upper.CL
##  Neutral - Pleasant          -2 0.4388537 27 -3.08810169 -0.9118983
##  Neutral - Unpleasant         1 0.4388537 27 -0.08810169  2.0881017
##  Pleasant - Unpleasant        3 0.4388537 27  1.91189831  4.0881017
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 3 estimates

Adjust2

##  contrast              estimate        SE df t.ratio p.value
##  Neutral - Pleasant          -2 0.4388537 27  -4.557  0.0002
##  Neutral - Unpleasant         1 0.4388537 27   2.279  0.0308
##  Pleasant - Unpleasant        3 0.4388537 27   6.836  <.0001
## 
## P value adjustment: holm method for 3 tests

Adjust3

##  contrast              estimate        SE df t.ratio p.value
##  Neutral - Pleasant          -2 0.4388537 27  -4.557  0.0003
##  Neutral - Unpleasant         1 0.4388537 27   2.279  0.0924
##  Pleasant - Unpleasant        3 0.4388537 27   6.836  <.0001
## 
## P value adjustment: bonferroni method for 3 tests

Adjust4

##  contrast              estimate        SE df t.ratio p.value
##  Neutral - Pleasant          -2 0.4388537 27  -4.557  0.0002
##  Neutral - Unpleasant         1 0.4388537 27   2.279  0.0308
##  Pleasant - Unpleasant        3 0.4388537 27   6.836  <.0001
## 
## P value adjustment: fdr method for 3 tests

La prueba de Bonferroni, muy versátil a la vez que sofocante, permitirá recordar la lógica a la base de este tipo de pruebas

---

Aclaraciones:
Tipos de p.adjust
“holm”, “hochberg”, “hommel”, “bonferroni”. Controlan “family-wise error rate”
“BH” o “fdr” (Benjamini & Hochberg). Control the false discovery rate (fdr)
“BY” (Benjamini & Yekutieli). También del tipo fdr.

Para obtener todos los tipos de contrastes que admite la librería
ls(“package:lsmeans”, pat=“.lsmc”)
pairwise.lsmc(levs, …)
revpairwise.lsmc(levs, …)
tukey.lsmc(levs, reverse = FALSE)
—
poly.lsmc(levs, max.degree = min(6, k - 1))
—
trt.vs.ctrl.lsmc(levs, ref = 1)
trt.vs.ctrl1.lsmc(levs, …)
trt.vs.ctrlk.lsmc(levs, …)
dunnett.lsmc(levs, ref = 1)
—
consec.lsmc(levs, reverse = FALSE, …)
mean_chg.lsmc(levs, reverse = FALSE, …)
—
eff.lsmc(levs, …)
del.eff.lsmc(levs, …)
—

Si el número de contraste no es elevado (e.g. para variables de hasta 3 ó 4 niveles) se recomendaría la prueba de Holm.
Por el contrario, si es muy elevado el nº de contrastes simultáneos (e.g. variables de 5 ó más nivles), se recomendaría Benjamini & Hochberg.
Un Ejemplo concreto:

Blanco, S., Hernández, R., Franchelli, G., Ramos-Álvarez, M.M., and Peinado, M.A. (2017). Melatonin influences NO/NOS pathway and reduces oxidative and nitrosative stress in a model of hypoxic-ischemic brain damage, Nitric Oxide, 62 (30), 32-43, http://dx.doi.org/10.1016/j.niox.2016.12.001.
Statistical significance was set at 0.05, and the probabilities obtained with the robust statistics were adjusted for multiple testing logic (pAdj onwards) for each of the sets of variables (clinical, oxidative, antioxidative, and health-related). Benjamini-Hochberg test was used, a less stringent condition than the usual Bonferroni variant, allowing gaining in statistical power (see Peña, Habiger, and Wu, 2011, for a review of the methods of multiple comparisons).

En este punto es muy recomendable la lectura del manual de Maxwell y Delaney (2004) con una nueva edición que saldrá en el 2017, o bien Ramos, Catena y Trujillo (2005) en castellano.
Para recapitular podría reinterpretar los resultados del Análisis Post Hoc a la luz de las nuevas pruebas

## Interpretación 3:  
##  
##  
##  
##  
##  
##  
##  
##  
##  
##  
###################################

4.5 Análisis para optimización del diseño

Estimaremos el Tamaño del Efecto y sus IC pero para los Contrastes del Diseño, pues a pesar de que cumpliríamos con las normas vigentes de de la APA (y lo que explícitamente incluyen buena parte de los comités de las revistas científicas) al incorporar estos cálculos únicamente en relación al ANOVA global, esto es como quedarse a medias.

# Requiere la librería {compute.es}
Grupo1=1; Grupo2=2 #Bastaría con ir cambiando esta parte
DatSel=Datos[Group %in% LosNivBetw[c(Grupo1:Grupo2)]]
LaT=t.test(DV~ Group,DatSel)$statistic[[1]]
n.1=c(DatSel[, .N, by = Group]$N[1])
n.2=c(DatSel[, .N, by = Group]$N[2])
#Para un contraste en particular
tes(t=LaT, n.1, n.2)

## Mean Differences ES: 
##  
##  d [ 95 %CI] = -2.12 [ -3.3 , -0.95 ] 
##   var(d) = 0.31 
##   p-value(d) = 0 
##   U3(d) = 1.69 % 
##   CLES(d) = 6.68 % 
##   Cliff's Delta = -0.87 
##  
##  g [ 95 %CI] = -2.03 [ -3.16 , -0.91 ] 
##   var(g) = 0.29 
##   p-value(g) = 0 
##   U3(g) = 2.11 % 
##   CLES(g) = 7.54 % 
##  
##  Correlation ES: 
##  
##  r [ 95 %CI] = 0.75 [ 0.42 , 0.9 ] 
##   var(r) = 0.01 
##   p-value(r) = 0 
##  
##  z [ 95 %CI] = 0.96 [ 0.45 , 1.47 ] 
##   var(z) = 0.06 
##   p-value(z) = 0 
##  
##  Odds Ratio ES: 
##  
##  OR [ 95 %CI] = 0.02 [ 0 , 0.18 ] 
##   p-value(OR) = 0 
##  
##  Log OR [ 95 %CI] = -3.85 [ -5.98 , -1.72 ] 
##   var(lOR) = 1.03 
##   p-value(Log OR) = 0 
##  
##  Other: 
##  
##  NNT = -5.04 
##  Total N = 20

# Si deseamos hacerlo de manera automatizada para todos los posibles contrastes
# Post Hoc, incluyendo mayor número de medidas (requiere las librerías {compute.es}
# y {WRS}), podrá valerse de la función propia CompESAPost (requiere cargar la
# librería {FuncMmRA.R}).
Result<-CompESAPost(Datos,"Group","DV")

## [1] "The Winsorized variance is zero for group 1"
## [1] "The Winsorized variance is zero for group 1"

Result

## [[1]]
##                 contr N.total n.1 n.2     d var.d  l.d   u.d U3.d cl.d
## 1 Neutral vs Pleasant      20  10  10 -2.12  0.31 -3.3 -0.95 1.69 6.68
##   cliffs.d pval.d     g var.g   l.g   u.g U3.g cl.g pval.g    r var.r  l.r
## 1    -0.87      0 -2.03  0.29 -3.16 -0.91 2.11 7.54      0 0.75  0.01 0.42
##   u.r pval.r fisher.z var.z  l.z  u.z   OR l.or u.or pval.or   lOR l.lor
## 1 0.9      0     0.96  0.06 0.45 1.47 0.02    0 0.18       0 -3.85 -5.98
##   u.lor pval.lor   NNT WilxoRob
## 1 -1.72        0 -5.04    0.785
## 
## [[2]]
##                   contr N.total n.1 n.2    d var.d  l.d  u.d  U3.d  cl.d
## 1 Neutral vs Unpleasant      20  10  10 1.13  0.23 0.12 2.15 87.16 78.87
##   cliffs.d pval.d    g var.g  l.g  u.g  U3.g  cl.g pval.g    r var.r  l.r
## 1     0.58   0.03 1.09  0.21 0.12 2.06 86.13 77.87   0.03 0.51  0.03 0.06
##    u.r pval.r fisher.z var.z  l.z  u.z   OR l.or  u.or pval.or  lOR l.lor
## 1 0.79   0.03     0.57  0.06 0.06 1.08 7.82 1.25 49.04    0.03 2.06  0.22
##   u.lor pval.lor  NNT WilxoRob
## 1  3.89     0.03 2.41    0.753
## 
## [[3]]
##                    contr N.total n.1 n.2    d var.d  l.d  u.d  U3.d  cl.d
## 1 Pleasant vs Unpleasant      20  10  10 2.71  0.38 1.41 4.02 99.67 97.25
##   cliffs.d pval.d   g var.g  l.g  u.g  U3.g  cl.g pval.g    r var.r  l.r
## 1     0.94      0 2.6  0.35 1.35 3.85 99.53 96.69      0 0.82  0.01 0.57
##    u.r pval.r fisher.z var.z  l.z  u.z     OR  l.or   u.or pval.or  lOR
## 1 0.93      0     1.16  0.06 0.65 1.66 137.27 12.94 1456.2       0 4.92
##   l.lor u.lor pval.lor NNT WilxoRob
## 1  2.56  7.28        0 1.3    0.855

#Para incluir la información Post Hoc en el apartado de resultados de nuestro
# informe (paper), al estilo clásico:
posthoc <- TukeyHSD(x=M.AOV)
Result$APA2<-ExtrTuk(posthoc) #requiere cargar la librería {FuncMmRA.R}
Result$APA2

## [[1]]
## [1] "Dif.Pleasant-Neutral=2; p= 3e-04; IC: [0.912, 3.088]"
## 
## [[2]]
## [1] "Dif.Unpleasant-Neutral=-1; p= 0.0762; IC: [-2.088, 0.088]"
## 
## [[3]]
## [1] "Dif.Unpleasant-Pleasant=-3; p= 0; IC: [-4.088, -1.912]"

• Pero sería preferible complementarlo con medidas del Tamaño Efecto como d ó g y sus correspondientes CI: Dif.Pleasant-Neutral=2; d=2.12; Vard= 0.31; 95% IC [0.95 3.30]; p= 3e-04. … Y así procedríamos con los otros 2 contrastes.
  
• Alternativamente, basado en el estadístico g (menor sesgo que d): Dif.Pleasant-Neutral=2; g=2.03; Vard= 0.29; 95% IC [0.91 3.16]; p= 3e-04. … Y así procedríamos con los otros 2 contrastes.
  
• Además podríamos añadir una interpretación de la magnitud neta de los estdísticos del Tamaño del Efecto (d o g), Mínimo= 0.41, Moderado = 1.15 y Fuerte = 2.70. Pero sería preferible poner en comparación su valor con el de otros estudios, en retrospectiva y/o en prospectiva.

---

# Finalmente, la estimación de la potencia (y así los cálculos del Tamaño del Efecto)
# tiene más sentido de manera prospectiva. Es decir, en el apartado del diseño
# podríamos estimar el Tamaño Muestral Óptimo para llegar a alcanzar un criterio prefijado
# de potencia estadística. 
# Recordemos los valores orientativos de la Tabla de Ferguson sobre la magnitud del efecto
# de Tratamiento, de manera que .04, .25 y .64, marcan respectivamente los valores mínimo,
# medio y elevado. Bastaría con incluir esta información en nuestro diseño particular.
pwr.anova.test(k=NNivBetw,f=.25,sig.level=.05,power=.8) #requiere la librería {pwr}

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 52.3966
##               f = 0.25
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: n is number in each group

En el ejemplo, el tamaño recomendado por cada uno de los 3 grupos de nuestro Diseño (recordar que se definieron en la variable “NNivBetw”) para lograr una potencia de 0.80 (valor muy razonable), manteniéndonos en el estándar de un 0.05 como nivel de significación, se estimaría en 53 sujetos aproximadamente.

En este punto es preferible obtener la familia de curvas de potencia variando los diferentes parámetros en los que se basa su cálculo, lo que haremos gracias a la función propia GraphPow desarrollada para este curso, que requiere tener cargada em memoria la librería {FuncMmRA.R},

GraphPow(NNivBetw,.05)

############################################## FIN ##############################################

Y aún mejor sería proyectar nuestras estimaciones de potencia en función de lo que otros investigadores han obtenido previamente, según la lógica de meta-análisis, para lo cual nos será de gran utilidad un gráfico tipo forest, bien a través de la función forest de R, o mejor a través de las librerías especializadas como {forestplot} o {metafor} …

---

5 Bases conceptuales y computacionales

---

Manuel Miguel Ramos Álvarez
Metodología de las Ciencias del Comportamiento
Dpto. Psicología, Universidad de Jaén, España
Grupo de investigación de Psicología Comparada
mramos@ujaen.es
http://www4.ujaen.es/~mramos

---