########################################################## ## Errores de tipo I y II, potencia de los tests, ## ## tipos de contrastes y sumas de cuadrados, ## ## problemas con interacciones que implican covariantes ## ########################################################## ## Luis M. Carrascal ## www.lmcarrascal.eu library(phia) ## para tests a posteriori y análisis de interacciones library(car) ## para obtener VIF's y usar Anova, bootCase, y transformación Box-Cox ###################################### ## para entender el ERROR DE TIPO I ## ###################################### ## muuuuuchos modelos nulos ## https://en.wikipedia.org/wiki/Type_I_and_type_II_errors ## eror de tipo I: RECHAZAR LA HIPÓTESIS NULA CUANDO DE HECHO ES CIERTA rm(list=ls()) ## con un caso concreto sencillo mirespuesta <- rnorm(n=80, mean=100, sd=25) factorA <- as.factor(c(rep(0, times=40), rep(1, times=40))) factorB <- as.factor(rep(c(0, 1), length.out=80)) datos <- data.frame(mirespuesta, factorA, factorB) modelo <- lm(mirespuesta ~ factorA*factorB, data=datos, contrasts=list(factorA=contr.sum, factorB=contr.sum)) summary(modelo) names(summary(modelo)) round(summary(modelo)$coefficients, 3) round(summary(modelo)$coefficients, 3)[2,4] round(summary(modelo)$coefficients, 3)[3,4] round(summary(modelo)$coefficients, 3)[4,4] ## interactionMeans(modelo) interactionMeans(modelo, factors="factorA") interactionMeans(modelo, factors="factorB") ## lo repetimos en un bucle con muchas repeticiones ## y además podemos modificar el tamaño muestral ## ¿qué tamaño de muestra quieres POR CELDA ## del diseño bifactorial con dos niveles por factor? mi_N <- 20 ## ¿qué alfa (P) deseas examinar? mi_P <- 0.05 pA <- rep(9999, times=10000) pB <- rep(9999, times=10000) pAB <- rep(9999, times=10000) for (i in 1:10000) { mirespuesta <- rnorm(n=mi_N*4, mean=100, sd=25) factorA <- as.factor(c(rep(0, times=mi_N*2), rep(1, times=mi_N*2))) factorB <- as.factor(rep(c(0, 1), length.out=mi_N*4)) datos <- data.frame(mirespuesta, factorA, factorB) modelo <- lm(mirespuesta ~ factorA*factorB, data=datos, contrasts=list(factorA=contr.sum, factorB=contr.sum)) pA[i] <- ifelse(summary(modelo)$coefficients[2,4]