library(multcomp)
library(nlme)
library(knitr)
load("~/Lola/GAB2020.RData")Los efectos fijos para cada uno de los parámetros debe compararse a partír de algún método que contemple controlar la tasa de error de tipo I. Para ello se propone el método de Benjamini- Hochberg. Se utilizó la libreria multcomp y se proporcionaron los contrastes necesarios para los cálculos a través de la matriz COM
Primero tomamos el vector de estimaciónes de efectos fijos
fixef(MODELO_101_d)Para luego crear por un lado un vector (M)que contiene las estimaciónes para el parámetro A para cada una de las temporadas
M <- rbind(c(1, 0, 0,rep(0,12)),
c(1, 1, 0,rep(0,12)),
c(1, 0, 1,rep(0,12)))
rownames(M) <- c("T1718", "T1819","T1920")
#MY un vector (K1) que contiene los contrastes entre las niveles de A para temporadas:
K1<-rbind(c(-1, 1, 0),
c(-1, 0, 1),
c(0, 1, -1))
rownames(K1) <- c(" A T1718 - T1819", " A T1718 - T1920"," A T1819 - T1920")Para obtener el producto de ambos vectores, el cual nos permite obtener el vector de contrastes a aplicar al vector de coeficientes del modelo:
COM<-K1%*%M
COMLa prueba para los niveles de A en TEMPORADAS finalmente se obtiene:
resulta<-glht(MODELO_101_d, linfct = COM)
summary(resulta,test=adjusted(type = "BH"))Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
A T1718 - T1819 == 0 -0.01718 0.01922 -0.894 0.371
A T1718 - T1920 == 0 -0.10295 0.01877 -5.485 1.24e-07 ***
A T1819 - T1920 == 0 0.08577 0.02024 4.239 3.37e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- BH method)De la misma forma se obtiene la prueba para los niveles de B en TEMPORADAS:
M <- rbind(c(rep(0,5),1, 0, 0,rep(0,7)),
c(rep(0,5),1, 1, 0,rep(0,7)),
c(rep(0,5),1, 0, 1,rep(0,7)))
rownames(M) <- c("T1718", "T1819","T1920")
#M
K1<-rbind(c(-1, 1, 0),
c(-1, 0, 1),
c(0, 1, -1))
rownames(K1) <- c("T1718 - T1819", " B T1718 - T1920"," B T1819 - T1920")
#K1
COM<-K1%*%M
resulta<-glht(MODELO_101_d, linfct = COM)
summary(resulta,test=adjusted(type = "BH"))Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
T1718 - T1819 == 0 -0.6510 0.1906 -3.416 0.000952 ***
B T1718 - T1920 == 0 -0.8148 0.1730 -4.710 7.42e-06 ***
B T1819 - T1920 == 0 0.1638 0.2016 0.813 0.416433
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- BH method)Y de C en TEMPORADAS:
M <- rbind(c(rep(0,10),1, 0, 0,rep(0,2)),
c(rep(0,10),1, 1, 0,rep(0,2)),
c(rep(0,10),1, 0, 1,rep(0,2)))
rownames(M) <- c("T1718", "T1819","T1920")
#M
K1<-rbind(c(-1, 1, 0),
c(-1, 0, 1),
c(0, 1, -1))
rownames(K1) <- c(" C T1718 - T1819", " C T1718 - T1920"," C T1819 - T1920")
#K1
COM<-K1%*%M
resulta<-glht(MODELO_101_d, linfct = COM)Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
C T1718 - T1819 == 0 -0.006753 0.004711 -1.433 0.151760
C T1718 - T1920 == 0 0.012816 0.004055 3.161 0.002362 **
C T1819 - T1920 == 0 -0.019568 0.005087 -3.847 0.000359 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- BH method)El mismo procedimiento se realiza para A en TAMAÑOS :
M <- rbind(c(1, 0, 0,0,0,rep(0,10)),
c(1, 0,0,1, 0,rep(0,10)),
c(1, 0,0,0, 1,rep(0,10)))
rownames(M) <- c("GRANDE", "MEDIANO","PEQUEÑO")
#M
K1<-rbind(c(-1, 1, 0),
c(-1, 0, 1),
c(0, 1, -1))
rownames(K1) <- c("A GRANDE-MEDIANO", "A GRANDE-PEQUEÑO","A MEDIANO-PEQUEÑO")
#K1
COM<-K1%*%M
resulta<-glht(MODELO_101_d, linfct = COM)
summary(resulta,test=adjusted(type = "BH"))Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
A GRANDE-MEDIANO == 0 0.08860 0.01914 4.630 3.67e-06 ***
A GRANDE-PEQUEÑO == 0 0.20905 0.01916 10.913 < 2e-16 ***
A MEDIANO-PEQUEÑO == 0 -0.12045 0.01914 -6.292 4.72e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- BH method)El mismo procedimiento se realiza para B en TAMAÑOS :
M <- rbind(c(rep(0,5),1, 0, 0,0,0,rep(0,5)),
c(rep(0,5),1, 0,0,1, 0,rep(0,5)),
c(rep(0,5),1, 0,0,0, 1,rep(0,5)))
rownames(M) <- c("T1718", "T1819","T1920")
#M
K1<-rbind(c(-1, 1, 0),
c(-1, 0, 1),
c(0, 1, -1))
rownames(K1) <- c("B GRANDE-MEDIANO", "B GRANDE-PEQUEÑO","B MEDIANO-PEQUEÑO")
#K1
COM<-K1%*%M
resulta<-glht(MODELO_101_d, linfct = COM)
summary(resulta,test=adjusted(type = "BH"))Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B GRANDE-MEDIANO == 0 0.9060 0.1766 5.132 2.87e-07 ***
B GRANDE-PEQUEÑO == 0 2.4799 0.1851 13.400 < 2e-16 ***
B MEDIANO-PEQUEÑO == 0 -1.5739 0.1882 -8.364 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- BH method)El mismo procedimiento se realiza para C en TAMAÑOS :
M <- rbind(c(rep(0,10),1, 0,0,0, 0),
c(rep(0,10),1, 0,0,1, 0),
c(rep(0,10),1, 0,0,0, 1))
rownames(M) <- c("T1718", "T1819","T1920")
#M
K1<-rbind(c(-1, 1, 0),
c(-1, 0, 1),
c(0, 1, -1))
rownames(K1) <- c("C GRANDE-MEDIANO", "C GRANDE-PEQUEÑO","C MEDIANO-PEQUEÑO")
#K1
COM<-K1%*%M
resulta<-glht(MODELO_101_d, linfct = COM)
summary(resulta,test=adjusted(type = "BH"))Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
C GRANDE-MEDIANO == 0 -0.013601 0.004291 -3.170 0.00152 **
C GRANDE-PEQUEÑO == 0 -0.040543 0.004341 -9.340 < 2e-16 ***
C MEDIANO-PEQUEÑO == 0 0.026942 0.004365 6.172 1.01e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- BH method)Benjamini Y., Hochberg Y.1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B, 57:289-300.
multcomp: Simultaneous Inference in General Parametric Models Versión 1.4-14. Paquete de R. Autores: Frank Brets, Peter Westfall. Colaboradores: Richard M. Heiberger, Andre Schuetzenmeister, Susan Scheibe