##################################### ## EL MISTERIO OCULTO TRAS LA P ## ## "SIGNIFICACIÓN" FRECUENTISTA ## ##################################### ## Luis M. Carrascal ## https://lmcarrascal.eu ## 05/05/2022 ## corrido bajo R-4.0.0 ## para ver qué tenemos y borrar datos rm(list=ls()) ## para limpiar la consola cat("\014") ## para limpiar todos los plots en el panel gráfico dev.off() ## para restablecer el generador aleatorio set.seed(NULL) ## cargamos librerías library(phia) library(MuMIn) library(MASS) library(psych) ## lecturas recomendadas ## https://www.aub.edu.lb/SHARP/PublishingImages/Pages/publications/proposal-lower-P-value-thresholds.pdf ## https://www.r-bloggers.com/whats-the-probability-that-a-significant-p-value-indicates-a-true-effect/ ## https://royalsocietypublishing.org/doi/full/10.1098/rsos.140216#xref-ref-2-1 ##################################### ## ERROR DE TIPO I CON CORRELACION ## ##################################### ## generamos los datos mi_N <- 15 variable1 <- rnorm(n=mi_N, mean=0, sd=1) variable2 <- rnorm(n=mi_N, mean=0, sd=1) ## hacemos la correlación aleatoria plot(variable1, variable2) cor(variable1, variable2) cor.test(variable1, variable2)$p.value ## lo repetimos 10000 veces mi_N <- 15 r.alazar_t1 <- rep(9999, times=10000) p.alazar_t1 <- rep(9999, times=10000) for (i in 1:10000) { variable1 <- rnorm(n=mi_N, mean=0, sd=1) variable2 <- rnorm(n=mi_N, mean=0, sd=1) r.alazar_t1[i] <- cor(variable1, variable2) p.alazar_t1[i] <- cor.test(variable1, variable2)$p.value } ## dev.off() plot(r.alazar_t1, p.alazar_t1, xlab="correlaciones aleatorias", ylab="significaciones (p) aleatorias") title(main="Error de tipo I", col.main="red", cex.main=2) hist(r.alazar_t1, breaks=15) abline(v=mean(r.alazar_t1), col="blue") abline(v=quantile(r.alazar_t1, probs=0.025), col="red") abline(v=quantile(r.alazar_t1, probs=0.975), col="red") hist(p.alazar_t1, breaks=15) percentile_p_t1 <- ecdf(p.alazar_t1) percentile_p_t1(0.05) ## error de tipo I (rechazar la Ho siendo correcta) 1-percentile_p_t1(0.05) ## aceptar la Ho siendo correcta mean(r.alazar_t1) ## effect size ###################################### ## ERROR DE TIPO II CON CORRELACION ## ###################################### ## generamos los datos asumiendo que sí existe correlación mi_N <- 15 variable1 <- rnorm(n=mi_N, mean=4, sd=1) variable2 <- 10 + 6*variable1 + rnorm(n=mi_N, mean=0, sd=9) ## hacemos la correlación que asume asociación plot(variable1, variable2) abline(a=10, b=6, col="red", lwd=2) cor(variable1, variable2) cor.test(variable1, variable2)$p.value ## lo repetimos 10000 veces mi_N <- 15 r.alazar_t2 <- rep(9999, times=10000) p.alazar_t2 <- rep(9999, times=10000) for (i in 1:10000) { variable1 <- rnorm(n=mi_N, mean=4, sd=1) variable2 <- 10 + 6*variable1 + rnorm(n=mi_N, mean=0, sd=9) r.alazar_t2[i] <- cor(variable1, variable2) p.alazar_t2[i] <- cor.test(variable1, variable2)$p.value } ## dev.off() plot(r.alazar_t2, p.alazar_t2, xlab="correlaciones aleatorias", ylab="significaciones (p) aleatorias") title(main="Error de tipo II", col.main="red", cex.main=2) hist(r.alazar_t2, breaks=15) abline(v=mean(r.alazar_t2), col="blue") abline(v=quantile(r.alazar_t2, probs=0.025), col="red") abline(v=quantile(r.alazar_t2, probs=0.975), col="red") hist(p.alazar_t2, breaks=15) percentile_p_t2 <- ecdf(p.alazar_t2) percentile_p_t2(0.05) ## potencia (1-beta) 1-percentile_p_t2(0.05) ## error de tipo II (beta) mean(r.alazar_t2) ## effect size ############################### ## ERROR DE TIPO I CON ANOVA ## ############################### ## generamos los datos tamañogrupo <- 15 factor01 <- as.factor(rep(c(0, 1), times=tamañogrupo)) normal.aleatoria <- rnorm(n=tamañogrupo*2, mean=25, sd=5) misdatos_t1 <- data.frame(factor01, normal.aleatoria) ## si tenemos un factor, correremos la siguiente línea para establecer unos contrastes correctos options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## hacemos el modelo aleatorio modeloaleatorio <- lm(normal.aleatoria ~ factor01, data=misdatos_t1) plot(interactionMeans(modeloaleatorio), legend.margin=0.2) summary(modeloaleatorio) ## summary(modeloaleatorio)$r.squared summary(modeloaleatorio)$coefficients summary(modeloaleatorio)$coefficients[2,4] ## hacemos el modelo aleatorio nulo modeloaleatorio.nulo <- update(modeloaleatorio, .~. -factor01) ## comparamos los dos modelos utilizando Akaike information criterion (AIC) AICc(modeloaleatorio.nulo, modeloaleatorio) deltaAICc <- AICc(modeloaleatorio) - AICc(modeloaleatorio.nulo) print(c("cuántas veces es mejor mimodelo =", exp(-0.5*deltaAICc)), quote=F) ## lo repetimos 10000 veces tamañogrupo <- 15 r2.alazar_t1 <- rep(9999, times=10000) p.alazar_t1 <- rep(9999, times=10000) delta.azar_t1 <- rep(9999, times=10000) for (i in 1:10000) { factor01 <- as.factor(rep(c(0, 1), times=tamañogrupo)) normal.aleatoria <- rnorm(n=tamañogrupo*2, mean=25, sd=5) misdatos_t1 <- data.frame(factor01, normal.aleatoria) modeloaleatorio <- lm(normal.aleatoria ~ factor01, data=misdatos_t1) modeloaleatorio.nulo <- update(modeloaleatorio, .~. -factor01) r2.alazar_t1[i] <- summary(modeloaleatorio)$r.squared p.alazar_t1[i] <- summary(modeloaleatorio)$coefficients[2,4] delta.azar_t1[i] <- AICc(modeloaleatorio) - AICc(modeloaleatorio.nulo) } ## dev.off() plot(r2.alazar_t1, p.alazar_t1) abline(h=0.05, col="red") abline(h=0.005, col="blue") ## plot(delta.azar_t1, p.alazar_t1) abline(h=0.05, col="red") abline(h=0.005, col="blue") ## plot(exp(-0.5*delta.azar_t1), p.alazar_t1) abline(h=0.05, col="red") abline(h=0.005, col="blue") ## hist(p.alazar_t1, breaks=15) percentile_p_t1 <- ecdf(p.alazar_t1) percentile_p_t1(0.05) ## error de tipo I 1-percentile_p_t1(0.05) ## aceptar la Ho siendo correcta ## hist(r2.alazar_t1, breaks=10) quantile(r2.alazar_t1, probs=0.95) ################################ ## ERROR DE TIPO II CON ANOVA ## ################################ factor01 <- as.factor(c(rep(0, times=5000), rep(1, times=5000))) respuesta <- c(rnorm(n=5000, mean=22, sd=5), rnorm(n=5000, mean=28, sd=5)) misdatos_t2 <- data.frame(factor01, respuesta) names(misdatos_t2) ## modelo omnisciente modelo.omnisciente <- lm(respuesta ~ factor01, data=misdatos_t2) summary(modelo.omnisciente) ## hacemos el modelo con una muestra concreta tamañogrupo <- 15 factor01 <- as.factor(c(rep(0, times=tamañogrupo), rep(1, times=tamañogrupo))) respuesta <- c(rnorm(n=tamañogrupo, mean=22, sd=5), rnorm(n=tamañogrupo, mean=28, sd=5)) datos_t2 <- data.frame(factor01, respuesta) names(datos_t2) modelo <- lm(respuesta ~ factor01, data=datos_t2) aggregate(respuesta ~ factor01, data=datos_t2, FUN=mean) plot(interactionMeans(modelo), legend.margin=0.2) summary(modelo) ## modelo nulo modelo.nulo <- update(modelo, .~. -factor01) AICc(modelo.nulo, modelo) deltaAICc <- AICc(modelo) - AICc(modelo.nulo) print(c("cuántas veces es mejor mimodelo =", exp(-0.5*deltaAICc)), quote=F) ## hacemos el experimento del error de tipo II 2000 veces tamañogrupo <- 15 r2.alazar_t2 <- rep(9999, times=10000) p.alazar_t2 <- rep(9999, times=10000) coeficiente_t2 <- rep(9999, times=10000) delta.azar_t2 <- rep(9999, times=10000) for (i in 1:10000) { factor01 <- as.factor(c(rep(0, times=tamañogrupo), rep(1, times=tamañogrupo))) respuesta <- c(rnorm(n=tamañogrupo, mean=22, sd=5), rnorm(n=tamañogrupo, mean=28, sd=5)) datos_t2 <- data.frame(factor01, respuesta) names(datos_t2) modelo <- lm(respuesta ~ factor01, data=datos_t2) modelo.nulo <- update(modelo, .~. -factor01) r2.alazar_t2[i] <- summary(modelo)$r.squared p.alazar_t2[i] <- summary(modelo)$coefficients[2,4] coeficiente_t2[i] <- summary(modelo)$coefficients[2,1] delta.azar_t2[i] <- AICc(modelo) - AICc(modelo.nulo) } ## dev.off() plot(r2.alazar_t2, p.alazar_t2) abline(h=0.05, col="red") abline(h=0.005, col="blue") ## plot(delta.azar_t2, p.alazar_t2) abline(h=0.05, col="red") abline(h=0.005, col="blue") ## plot(exp(-0.5*delta.azar_t2), p.alazar_t2) abline(h=0.05, col="red") abline(h=0.005, col="blue") ## hist(p.alazar_t2, breaks=15) percentile_p_t2 <- ecdf(p.alazar_t2) percentile_p_t2(0.05) ## potencia (1-beta) 1-percentile_p_t2(0.05) ## error de tipo II (beta) ## hist(coeficiente_t2, breaks=10) percentile_coef_t2 <- ecdf(coeficiente_t2) percentile_coef_t2(0) ## hist(r2.alazar_t2, breaks=10) quantile(r2.alazar_t2, probs=0.05) ########################################## ## ERRORES DE TIPO I Y II CON REGRESION ## ########################################## set.seed(1) x_e2 <- rnorm(n=5000, mean=100, sd=20) ## predictora con mucha varianza asociada al error de tipo II z_e2 <- rnorm(n=5000, mean=100, sd=8) ## predictora con poca varianza asociada al error de tipo II v_e2 <- rnorm(n=5000, mean=100, sd=20) ## predictora con mucha varianza asociada al error de tipo II w_e1 <- rnorm(n=5000, mean=100, sd=20) ## predictora con mucha varianza asociada al error de tipo I ruido <- rnorm(n=5000, mean=0, sd=20) ## ruido introducido al generar la variable respuesta y ## generamos la variable respuesta y ## x_e2 y z_e2 tienen el mismo coeficiente (magnitud de efecto) ## v_e2 tiene menor coeficiente que x_e2 y z_e2 (magnitud de efecto menor) ## w_e1 establecemos que no se relaciona con y al asignarle un coeficiente 0 y <- 100 + 1*x_e2 + 1*z_e2 + 0.5*v_e2 + 0*w_e1 + ruido ## creamos la matriz de datos (dataframe) con las variables "sueltas" generadas previamente ## lo denominamos "población" de estudio, que posteriormente vamos a "muestrear" poblacion <- data.frame(y, x_e2, z_e2, v_e2, w_e1) ## representamos las relaciones entre ellas pairs(~y+x_e2+v_e2+z_e2+w_e1, data=poblacion) ## estimamos el modelo de relaciones entre respuesta y predictores ## esta es la realidad que pretendemos estudiar omnisciente <- lm(y ~x_e2+v_e2+z_e2+w_e1, data=poblacion) summary(omnisciente) ## vamos a muestrear la población original extrayendo unas observaciones (filas) al azar ## id contiene una selección al azar de esas filas ## en el argumento size establecemos nuestra muestra (sample-size) id <- sample(1:5000, replace=F, size=20) ## obtenemos la matriz de datos con nuestra muestra ## resultado de la actividad investigadora muestra <- poblacion[id, ] ## analizamos nuestra muestra y vemos cómo sus resultados coinciden o no con ## los "omnisicientes" que se dan en toda la población modelo <- lm(y ~x_e2+v_e2+z_e2+w_e1, data=muestra) ## resultados exhaustivos summary(modelo) ## versión abreviada de la tabla de coeficientes ## nos fijamos si la P es <0.05 round(summary(modelo)$coefficients, 3) ## repetimos el proceso varias veces corriendo estas líneas ## valorando como aparecen significaciones falsas (error de tipo I) con +w_e1 ## valorando que no aparecen algunas veces resultados significativos para x_e2, v_e2 y z_e2 (errores de tipo II) id <- sample(1:5000, replace=F, size=20) muestra <- poblacion[id, ] modelo <- lm(y ~x_e2+v_e2+z_e2+w_e1, data=muestra) round(summary(modelo)$coefficients, 3)