############################## ## MIS PRIMEROS PASOS CON R ## ############################## ## Luis M. Carrascal ## www.lmcarrascal.eu ## 24 Enero de 2017 ################################ ## DISTRIBUCIONES Y VARIABLES ## ################################ ## para generar una variable al azar según una distribución CONTINUA determinada var.uniforme <- runif(n=100, min=1, max=4) ## distribución uniforme var.normal <- rnorm(n=100, mean=2.5, sd=1) ## distribución gausiana var.poisson <- rpois(n=100, lambda=2.5) ## distribución Poisson var.negbin <- rnbinom(n=100, size=1, mu=2.5) ## distribución Binomial Negativa ## para generar una variable al azar según una distribución NOMINAL determinada var.binomial2 <- rbinom(n=100, prob=0.5, size=1) ## distrib. binomial con 0-1 var.binomial3 <- rbinom(n=100, prob=0.5, size=2) ## distrib. binomial con 0-2 length(var.normal) mean(var.normal) sd(var.normal) varianza.normal <- sd(var.normal)^2 varianza.normal ## identidad de las variables class(var.uniforme) class(var.poisson) class(var.binomial3) class(varianza.normal) ls() rm(varianza.normal) ls() ## histogramas sencillos hist(var.normal) hist(var.binomial2) hist(var.binomial3) ## gráficas y comparación de distribuciones par(mfcol=c(1,1)) par(mfrow=c(2,2)) ## fija los paneles según 2 columna y 2 filas hist(var.uniforme, density=5, freq=FALSE, main="distribución UNIFORME", xlim=c(0,10), ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal), sd=sd(var.normal)), col="red", lwd=2, add=TRUE, yaxt="n") hist(var.normal, density=5, freq=FALSE, main="distribución NORMAL", xlim=c(0,10), ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal), sd=sd(var.normal)), col="red", lwd=2, add=TRUE, yaxt="n") hist(var.poisson, density=5, freq=FALSE, main="distribución POISSON", xlim=c(0,10), ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal), sd=sd(var.normal)), col="red", lwd=2, add=TRUE, yaxt="n") hist(var.negbin, density=5, freq=FALSE, main="distribución BINOMIAL NEGATIVA", xlim=c(0,10), ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal), sd=sd(var.normal)), col="red", lwd=2, add=TRUE, yaxt="n") par(mfcol=c(1,1)) par(mfcol=c(1,1)) par(mfcol=c(1,2)) ## fija los paneles según 1 columna y 2 filas hist(var.normal, density=5, freq=FALSE, main="distribución simulada NORMAL", ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal), sd=sd(var.normal)), col="red", lwd=2, add=TRUE, yaxt="n") qqnorm(var.normal) qqline(var.normal, col="red", lwd=2) par(mfcol=c(1,1)) ## con un bucle para n=50 for (i in 1:10) { var.normal.50 <- rnorm(n=50, mean=2.5, sd=1) par(mfcol=c(1,1)) par(mfcol=c(1,2)) hist(var.normal.50, density=5, freq=FALSE, main="distribución simulada", ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal.50), sd=sd(var.normal.50)), col="red", lwd=2, add=TRUE, yaxt="n") qqnorm(var.normal.50) qqline(var.normal.50, col="red", lwd=2) par(mfcol=c(1,1)) } ## convertir un vector numérico con decimales a un vector de números enteros ## trampa; no lo hagais para sindecimales <- as.integer(var.uniforme) class(sindecimales) class(var.uniforme) ## para definir un factor a partir de un vector fbinom2 <- as.factor(var.binomial2) class(fbinom2) class(var.binomial2) ## para convertir una variable continua a otra con valores (0 1) ## a la condición "<2.5" se le asigna el valor 0, y a lo demás el valor 1 mi.factor2 <- ifelse(var.uniforme<2.5, 0, 1) ## y ahora la convertimos en factor mi.factor2 <- as.factor(mi.factor2) ## asignar valores de un vector a códigos-niveles de un factor ## considerando los puntos de corte incluidos en c(. . .) mi.factor3 <- as.factor(cut(var.uniforme, c(-1,2,3,5), labels=c("POCO", "MEDIO", "MUCHO"))) class(mi.factor3) ## es lo mismo porque las etiquetas son textuales y no pueden ser números mi.factor3 <- cut(var.uniforme, c(-1,2,3,5), labels=c("POCO", "MEDIO", "MUCHO")) class(mi.factor3) ## generamos variables manualmente mi.variable <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15) mi.variable <- c(1:15) mi.variable.seq <- seq(from=1, to=15, by=1) otra.variable <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) otra.variable.rep <- rep(0, times=20) otra.variable.rep <- rep(0, length.out=20) otra.variable.rep2 <- rep(c(1, 2, 3), times=5) ## combinar dos vectores en uno, añadiendo el segundo a continuación del primero mi.primera <- rep(0, times=20) mi.segunda <- rep(1, times=20) mi.primera.segunda <- c(mi.primera, mi.segunda) mi.segunda.primera <- c(mi.segunda, mi.primera) ## ordeno una variable "suelta" var.normal.ordenada <- sort(var.normal) ## ¡qué horror! plot(var.normal, var.normal.ordenada) plot(var.normal, sort(var.normal)) ## hacemos nuestra primera matriz de datos ## deben tener la misma longitud TODOS los vectores datos <- data.frame(var.uniforme, var.normal, var.poisson, var.negbin, var.binomial2, var.binomial3) ## nombres de las columnas ó variables names(datos) colnames(datos) ## clase de objeto, tipo de datos y dimensiones class(datos) str(datos) dim(datos) dim(datos)[1] ## número de filas dim(datos)[2] ## número de columnas ó variables ## ordeno mi matriz de datos (data.frame) por las filas de una variable, creando otra nueva matriz datos.ordenados <- datos[order(datos$var.uniforme),] ## cambio los nombres de las variables names(datos) <- c("uniforme", "normal", "poisson", "binomial_negativa", "binomial_01", "trinomial_0123") names(datos) ## para seleccionar-extraer una parte de mis casos-observaciones; dos formas equivalentes datos2 <- datos[c(1:15, 22:30),] datos3 <- subset(datos[c(1:15, 22:30),]) ## más complejo: igual; mayor o igual O (|, AltGr+1) menor o igual; que cumplan las dos condiciones unidas por & datos4 <- subset(datos, binomial_01 == 0) datos5 <- subset(datos, normal >= 3 | normal <= 1.5) datos6 <- subset(datos, normal > 2 & poisson <= 2) #################################### ## MIS PRIMEROS PASOS CON MODELOS ## #################################### ## generamos variables asociadas para un modelo inventado ## para obtener SIEMPRE los mismos valores simulados fijamos un valor "seed" set.seed(2009) x <- rnorm(n=50, mean=30, sd=10) z <- rnorm(n=50, mean=30, sd=10) ruido <- rnorm(n=50, mean=0, sd=50) y <- 200 + 5*x - 5*z + ruido ## obtenemos dos gráficas (id de una a otra con la flecha izda y dcha de la esquina superior del panel PLOTS) hist(y); qqnorm(y) ## y restablecemos la aleatoridad no controlada por nosotros set.seed(NULL) ## al correrlo de este modo sin set.seed(...), siempre se obtienen resultados diferentes ## seleccionamos el bloque y ctrl+intro ## comparamos como cambian las cosas; antes no x <- rnorm(n=50, mean=30, sd=10) z <- rnorm(n=50, mean=30, sd=10) ruido <- rnorm(n=50, mean=0, sd=50) y <- 200 + 5*x - 5*z + ruido qqnorm(y) ## hacemos un factor var.binom <- rbinom(n=50, prob=0.5, size=1) factor <- as.factor(cut(var.binom, c(-1,0,1), labels=c("CONTROL", "TRATAMIENTO"))) ## creamos un data.frame a partir de variables "sueltas" datos <- data.frame(y, x, z, ruido, factor) ## obtenemos unos valores descriptivos library(psych) describe(datos) mi.tabla <- as.data.frame(describe(datos)) mi.tabla mi.tabla.reducida <- mi.tabla[,c(2:5, 8:9, 10:12)] mi.tabla.reducida round(mi.tabla.reducida, 3) ## hacemos una figura sencilla, pero elaborándola bastante ## http://www.statmethods.net/advgraphs/parameters.html ## pch es el tipo de punto; cex es su tamaño plot(datos$x, datos$y, pch=1, cex=1.5, xlab="predictor x", ylab="respuesta y", col.lab="red", main="MI GRÁFICO DE RELACIÓN X - Y", col.main="purple") abline(v=mean(datos$x), col="green") ## media de la x abline(h=mean(datos$y), col="green") ## media de la y abline(lm(datos$y~datos$x), col="blue", lty=5, lwd=3) ## hacemos una figura compleja pero sin elaborar eqt1 <- as.formula(y ~ x + z + ruido) eqt2 <- as.formula(y ~ x + z) pairs(eqt1, data=datos, main="mi figura de relaciones entre variables") ## hacemos un modelo NULO modelo.nulo <- lm(z ~ x + factor, data=datos) names(modelo.nulo) modelo.nulo$model summary(modelo.nulo) names(summary(modelo.nulo)) summary(modelo.nulo)$r.squared summary(modelo.nulo)$coefficients summary(modelo.nulo)$coefficients[2,4] ## hacemos un modelo modelo <- lm(y ~ x + z, data=datos) names(modelo) modelo$model summary(modelo) names(summary(modelo)) summary(modelo)$r.squared summary(modelo)$coefficients summary(modelo)$coefficients[2,4] ## para entender cómo funcionan los bucles ## https://www.datacamp.com/community/tutorials/tutorial-on-loops-in-r#gs.On0i=QE ############################################## ## PARA ENTENDER LOS ERRORES DE TIPO I Y II ## ############################################## ## https://en.wikipedia.org/wiki/Type_I_and_type_II_errors ## ERROR DE TIPO I ## muuuuuchos modelos nulos px <- seq(1, 10000, 1) pfactor <- seq(1, 10000, 1) r2.nulo <- rep(0, times=10000) ## for (i in 1:10000) { x <- rnorm(n=50, mean=30, sd=10) z <- rnorm(n=50, mean=30, sd=10) var.binom <- rbinom(n=50, prob=0.5, size=1) factor <- as.factor(cut(var.binom, c(-1,0,1), labels=c("CONTROL", "TRATAMIENTO"))) modelo.nulo <- lm(z ~ x + factor) px[i] <- summary(modelo.nulo)$coefficients[2,4] pfactor[i] <- summary(modelo.nulo)$coefficients[3,4] r2.nulo[i] <- summary(modelo.nulo)$r.squared } np05.x <- ifelse(px<0.05, 1, 0) np05.pfactor <- ifelse(pfactor<0.05, 1, 0) prop05.x <- sum(np05.x)/10000 prop05.pfactor <- sum(np05.pfactor)/10000 r2.nulo.media <- mean(r2.nulo) ## los resultados de 10000 procesos son r2.nulo.media*100 ## media de los R2 de los modelos en % prop05.x ## error de tipo I para la predictora x prop05.pfactor ## error de tipo I para la predictora nominal factor ## ERROR DE TIPO II ## muuuuuchos modelos ALTERNATIVOS ## ... te pregunto: ¿cuánto ruido quieres? y ¿qué tamaño muestral? cantidad.de.ruido <- 75 mi.muestra <- 30 px <- seq(1, 10000, 1) pz <- seq(1, 10000, 1) r2.modelo <- rep(0, times=10000) ## for (i in 1:10000) { x <- rnorm(n=mi.muestra, mean=30, sd=10) z <- rnorm(n=mi.muestra, mean=30, sd=10) ruido <- rnorm(n=mi.muestra, mean=0, sd=cantidad.de.ruido) y <- 200 + 5*x - 5*z + ruido modelo <- lm(y ~ x + z) px[i] <- summary(modelo)$coefficients[2,4] pz[i] <- summary(modelo)$coefficients[3,4] r2.modelo[i] <- summary(modelo)$r.squared } np.no05.x <- ifelse(px<0.05, 0, 1) np.no05.z <- ifelse(pz<0.05, 0, 1) prop.no05.x <- sum(np.no05.x)/10000 prop.no05.z <- sum(np.no05.z)/10000 r2.modelo.media <- mean(r2.modelo) ## los resultados de 10000 procesos son r2.modelo.media*100 ## media de los R2 de los modelos en % prop.no05.x ## error de tipo II para la predictora x prop.no05.z ## error de tipo II para la predictora z ################################################### ## BOOTSTRAPPING DE UN MODELO GENERAL LIN. MODEL ## ################################################### ## http://www.uvm.edu/~dhowell/StatPages/Permutation%20Anova/PermTestsAnova.html ## limpiamos completamente el "Environment" para cargar un nuevo juego de datos y hacer otros análisis rm(list=ls()) ## obtener los datos de internet y ponderlos en uso en "environment" meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayMANOVA.RData") ## abro una conexión con internet load(meconecto) ## cargo el archivo de la conexión close(meconecto); rm(meconecto) ## cierro la conexión library(car) library(psych) options(contrasts=c(factor="contr.sum", ordered="contr.poly")) eqt <- as.formula(MUSCULO ~ TARSO + ZONA*SEXO) modelo <- lm(eqt, data=datos) nterms <- summary(modelo)$df[1] F.boot <- matrix(nrow=1000, ncol=nterms) p.boot <- matrix(nrow=1000, ncol=nterms) colnames(F.boot) <- rownames(Anova(modelo, type=3, test="F"))[c(1:nterms)] colnames(p.boot) <- rownames(Anova(modelo, type=3, test="F"))[c(1:nterms)] set.seed(111) for (i in 1:1000) { cat("\014"); print(i) ## para ver en qué proceso estoy ind <- sample(nrow(datos), size=nrow(datos), replace=TRUE) boot.datos <- datos[ind,] modelo.boot <- lm(eqt, data=boot.datos) F.boot[i,] <- Anova(modelo.boot, type=3, test="F")[c(1:nterms),3] p.boot[i,] <- Anova(modelo.boot, type=3, test="F")[c(1:nterms),4] } set.seed(NULL) F.boot <- as.data.frame(F.boot) p.boot <- as.data.frame(p.boot) ## estimas de significación de la F vs Ho de la F ## discusión sobre la hipótesis nula para la distribución F ## https://www.statlect.com/probability-distributions/F-distribution ## http://stats.stackexchange.com/questions/114174/in-anova-if-the-null-hypothesis-is-true-is-the-expected-value-of-f-guaranteed-t 1-sum(ifelse(F.boot[,2]<(modelo$df.residual/(modelo$df.residual-2)), 0, 1)/1000) ecdf(F.boot[,2])(modelo$df.residual/(modelo$df.residual-2)) ## calculamos las nuevas significaciones tabla.F <- matrix(nrow=nterms, ncol=4) colnames(tabla.F) <- c("F_modelo", "p_modelo", "p_nueva", "%_NOsignif_0.05") tabla.F <- as.data.frame(tabla.F) tabla.F[,1] <- Anova(modelo, type=3, test="F")[c(1:nterms),3] tabla.F[,2] <- Anova(modelo, type=3, test="F")[c(1:nterms),4] for (j in 1:nterms) { tabla.F[j,3] <- ecdf(F.boot[,j])(modelo$df.residual/(modelo$df.residual-2)) tabla.F[j,4] <- 100*(sum(ifelse(p.boot[,j]<0.05, 0, 1))/1000) hist(p.boot[,j], breaks=20, ylim=c(0, 1000), xlab="significación (p)", main=names(p.boot)[j]) } rownames(tabla.F) <- rownames(Anova(modelo, type=3, test="F"))[c(1:nterms)] Anova(modelo, type=3, test="F") round(tabla.F, 4) qf(0.05, 1, 33, lower.tail=FALSE) pf(1, 1, 33, lower.tail=FALSE) lapply(F.boot, quantile, probs=c(0.9, 0.95, 0.99, 0.999), name=TRUE) ## valores nuevos de significación recalculados a partir de F obtenidas de modelos nulos ## generamos 2.000 modelos nulos n.efectos <- length(modelo$coefficients) nef <- n.efectos F.nulas <- matrix(nrow=2000, ncol=nef) for (i in 1:2000) { analizar.nulos <- modelo$model analizar.nulos[,1] <- sample(modelo$model[,1]) modelo.nulo <- lm(eqt, data=analizar.nulos, contrasts=c(factor="contr.sum", ordered="contr.poly")) F.esta.vez <- Anova(modelo.nulo, type=3, test="F")[c(1:nef),3] F.nulas[i,] <- F.esta.vez } F.nulas <- as.data.frame(F.nulas) colnames(F.nulas) <- rownames(Anova(modelo.nulo, type=3, test="F"))[c(1:nef)] ## calculamos los cuantiles de las F nulas tabla.quantiles <- matrix(nrow=nef, ncol=5) colnames(tabla.quantiles) <- c("Fp_0.1", "Fp_0.05", "Fp_0.01", "Fp_0.005", "Fp_0.001") for (j in 1:nef) { tabla.quantiles[j,] <- quantile(F.nulas[,j], probs=c(0.9, 0.95, 0.99, 0.995, 0.999), na.rm=TRUE) } tabla.quantiles <- as.data.frame(tabla.quantiles) tabla.quantiles$F_observada <- Anova(modelo, type=3, test="F")[c(1:nef),3] rownames(tabla.quantiles) <- rownames(Anova(modelo, type=3, test="F"))[c(1:nef)] round(tabla.quantiles, 3) ## calculamos las nuevas significaciones tabla.F <- matrix(nrow=nef, ncol=3) colnames(tabla.F) <- c("F_modelo", "p_modelo", "p_nueva") tabla.F <- as.data.frame(tabla.F) tabla.F[,1] <- Anova(modelo, type=3, test="F")[c(1:nef),3] tabla.F[,2] <- Anova(modelo, type=3, test="F")[c(1:nef),4] for (j in 1:nef) { tabla.F[j,3] <- 1-ecdf(F.nulas[,j])(tabla.F[j,1]) } rownames(tabla.F) <- rownames(Anova(modelo, type=3, test="F"))[c(1:nef)] round(tabla.F, 5) ################# ## OTRAS COSAS ## ################# ## una "NO-NORMAL" que no nos importa que no lo sea ... ¡¡Y LO HAREMOS BIEN!! v1 <- rnorm(20, 5, 2) v2 <- rnorm(20, 20, 2) s1 <- rep(0, 20) s2 <- rep(1, 20) respuesta <- c(v1, v2) factor <- c(s1, s2) class(factor) factor <- as.factor(factor) mmm <- lm(respuesta ~ factor) qqnorm(respuesta, datax=TRUE) qqline(respuesta, datax=TRUE) qqnorm(mmm$residuals, datax=TRUE) qqline(mmm$residuals, datax=TRUE) ## EFECTO DEL TIPO DE CONTRASTES Y DE SS ## de internet y ponderlo en uso en "environment" meconecto <- url("http://www.lmcarrascal.eu/cursos/nwayMANOVA.RData") ## abro una conexión con internet load(meconecto) ## cargo el archivo de la conexión close(meconecto); rm(meconecto) ## cierro la conexión library(car) library(phia) modelo.tre <- lm(MUSCULO~TARSO+ZONA*SEXO*EDAD, data=datos) modelo.sum <- lm(MUSCULO~TARSO+ZONA*SEXO*EDAD, data=datos, contrasts=list(ZONA=contr.sum, SEXO=contr.sum, EDAD=contr.sum)) modelo.sum2 <- lm(MUSCULO~EDAD*ZONA*SEXO+TARSO, data=datos, contrasts=list(ZONA=contr.sum, SEXO=contr.sum, EDAD=contr.sum)) Anova(modelo.tre, type=3, test="F") Anova(modelo.sum, type=3, test="F") summary(modelo.tre) summary(modelo.sum) plot(interactionMeans(modelo.tre), legend.margin=0.2) plot(interactionMeans(modelo.sum), legend.margin=0.2) Anova(modelo.sum, type=3, test="F") Anova(modelo.sum2, type=3, test="F") anova(modelo.sum) anova(modelo.sum2) contr.des <- matrix(c(3,-1,-1,-1, -1,3,-1,-1, -1,-1,3,-1), ncol=3) datos$DOMINANCIA2 <- ifelse(datos$DOMINANCIA<3, 0, 1) datos$FDOMINANCIA <- as.ordered(datos$DOMINANCIA) datos$FDOMINANCIA2 <- cut(datos$DOMINANCIA, c(0,1,2,3,4), labels=c("MUY SUBORDINADO", "SUBORDINADO", "DOMINANTE", "MUY DOMINANTE")) modelo2.sum <- lm(MUSCULO~TARSO+FDOMINANCIA2, data=datos, contrasts=list(FDOMINANCIA2=contr.sum)) modelo2.des <- lm(MUSCULO~TARSO+FDOMINANCIA2, data=datos, contrasts=list(FDOMINANCIA2=contr.des)) summary(modelo2.sum) summary(modelo2.des) contr.sum(4) contr.des Anova(modelo2.sum, type=3, test="F") Anova(modelo2.des, type=3, test="F") plot(interactionMeans(modelo2.sum), legend.margin=0.2) contr.lin <- matrix(c(-3,-1,1,3), ncol=1) modelo2.lin <- lm(MUSCULO~TARSO+FDOMINANCIA2, data=datos, contrasts=list(FDOMINANCIA2=contr.lin)) summary(modelo2.lin) summary(modelo2.des) Anova(modelo2.lin, type=3, test="F") Anova(modelo2.des, type=3, test="F") contr.lin2 <- matrix(c(-3,-3,1,5), ncol=1) modelo2.lin2 <- lm(MUSCULO~TARSO+FDOMINANCIA2, data=datos, contrasts=list(FDOMINANCIA2=contr.lin2)) summary(modelo2.lin) summary(modelo2.lin2) Anova(modelo2.lin, type=3, test="F") Anova(modelo2.lin2, type=3, test="F") ## Enoooorme problema con las interacciones QUE IMPLICAN COVARIANTES x <- rnorm(n=30, mean=30, sd=10) z <- rnorm(n=30, mean=30, sd=10) ruido <- rnorm(n=30, mean=0, sd=75) y <- 200 + 5*x - 5*z + ruido modelo <- lm(y ~ x + z) pairs(y ~ x + z + ruido) ## AÑADIMOS LA INTERACCIÓN (de las covariantes x - z con x:z) modelo2 <- lm(y ~ x + z + x:z) ## y vemos cómo cambian las significaciones summary(modelo) summary(modelo2)