#################################### ## Modelos Generalizados Lineales ## ## MODELOS DE AUSENCIA DE EFECTOS ## ## RANDOMIZATION TESTS ## #################################### ## Luis M. Carrascal ## https://lmcarrascal.eu library(psych) ## elegimos un Modelo Generalizado Lineal modelo <- modelo.nb ## los datos del dataframe data (=datos) con los que hemos trabajado están en modelo$model ## la primera variable se corresponde con la respuesta, y las siguientes son los predictores ## extraemos en un vector la respuesta de nuestro modelo (podemos verla en eqt) respuesta <- modelo$model[,1] ## el tamaño muestral (número de observaciones) podemos verlo con dim(modelo$model) ## observaciones variables ## podemos crear un vector indice de observaciones, numerandolas de 1 hasta ... dim(modelo$model)[1] indice <- c(1:dim(modelo$model)[1]) ## esta secuencia numerica ordenada podemos aleatorizarla indice.azar <- sample(indice, replace=FALSE) ## utilizando ese indice aleatorizado podemos generar la respuesta "desordenada" (aleatorizada en la posición de los casos sin seguir el orden inicial) respuesta.azar <- respuesta[indice.azar] ## vemos que la respuesta original y la respuesta "barajada al azar" no se relacionan plot(respuesta, respuesta.azar) cor(respuesta, respuesta.azar) ## ahora contruimos un modelo aleatorio, desacoplando los valores de la respuesta de los que toman las predictoras modelo.azar <- update(modelo, respuesta.azar ~.) ## en esta tabla se obtienen los coeficientes y sus significaciones coeftest(modelo.azar, df=modelo.azar$df.residual) ## lo comparamos con el modelo real coeftest(modelo, df=modelo$df.residual) ## RANDOMIZATION TESTS o PERMUTATION TESTS ## https://en.m.wikipedia.org/wiki/Resampling_(statistics) ## https://www.uvm.edu/~statdhtx/StatPages/Randomization%20Tests/RandomizationTestsOverview.html ## https://www.rdatagen.net/post/permutation-test-for-a-covid-19-pilot-nursing-home-study/ ## Este proceso podemos repetirlo muchas veces para obtener coeficientes que encierran efectos nulos ## al desacoplar los valores de la respuesta de los de las variables predictoras. ## Podemos utilizar este procedimiento para construir nuestras significaciones en el caso de que nuestros modelos de interés ## no se hayan ajustado suficientemente bien a los supuestos canonicos de los modelos GLM ## (normalidad de los residuos, heterocedasticidad, etc) ## de manera que no utilizaremos las significaciones derivadas de las Z, Chi2, F, t ## Simulaciones de modelos que encierran efectos nulos n.simulaciones <- 2000 n.coeficientes <- length(modelo$coefficients) coeficientes.nulos <- as.data.frame(matrix(999999, nrow=n.simulaciones, ncol=n.coeficientes)) colnames(coeficientes.nulos) <- names(modelo.azar$coefficients) for (i in 1:n.simulaciones) { respuesta <- modelo$model[,1] indice.azar <- sample(c(1:dim(modelo$model)[1]), replace=FALSE) respuesta.azar <- respuesta[indice.azar] modelo.azar <- update(modelo, respuesta.azar ~.) coeficientes.nulos[i,] <- modelo.azar$coefficients } ## obtenemos los resultados de miles de coeficientes nulos ## la desviación típica (sd) en este tipo de procesos se corresponde con el error estándard (se) de los parámetros estimados ## se proporcionan los cuantiles a alfa=0.05 (95%, Q0.025 Q0.975) y alfa=0.01 (99%, Q0.005 Q0.995) tabla.coefs.nulos <- describe(coeficientes.nulos, quant=c(0.025, 0.975, 0.005, 0.995))[,c(2:4,11:12, 14:17)] ## añadimos los coeficientes reales del modelo de interés tabla.coefs.nulos$coeficientes_observados <- modelo$coefficients tabla.coefs.nulos$Z <- (tabla.coefs.nulos$coeficientes_observados - tabla.coefs.nulos$mean) / tabla.coefs.nulos$sd significacion <- (1-pt(abs(tabla.coefs.nulos$Z), df=999999, lower.tail = TRUE))*2 tabla.coefs.nulos$P <- significacion ## para que los coeficientes reales observados sean significativos deberían quedar fuera de los intervalos de confianza a esas alfas print(tabla.coefs.nulos[, c(1:3, 6:10)], digits=5) ## podríamos utilizar las P derivadas del estadístico Z si los valores de sesgo y kurtosis estan proximos a cero print(tabla.coefs.nulos[, c(4:5,11:12)], digits=5) ## tambien podemos visualizar los resultados nulos y los coeficientes observados ## si la linea verde que denota los coeficientes observados queda dentro de las lineas azules ## que definen los intervalos de confianza al 95% de los coeficientes nulos, ## entonces ese coeficiente observado no es significativo a alfa=0.05 for (i in 1:n.coeficientes) { hist(coeficientes.nulos[,i], main=colnames(coeficientes.nulos)[i], col.main="red", cex.main=2) abline(v=modelo$coefficients[i], col="green", lwd=3) abline(v=quantile(coeficientes.nulos[,i], probs=0.025), col="blue") abline(v=quantile(coeficientes.nulos[,i], probs=0.975), col="blue") }