############################## ## MIS PRIMEROS PASOS CON R ## ############################## ## Luis M. Carrascal ## www.lmcarrascal.eu ## Febrero 2019 ################################################# ## APRENDIENDO A MANEJAR VECTORES Y DATAFRAMES ## ################################################# ## 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, 4), 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 <- rnorm(n=25, mean=50, sd=20) var.normal.ordenada <- sort(var.normal) ## ¡qué horror! plot(var.normal, var.normal.ordenada) ## para ver qué tenemos y borrar datos ls() rm(var.normal.ordenada) ls() rm(list=ls()) ## para limpiar la consola cat("\014") ## para limpiar todos los plots en el panel gráfico dev.off() ################################ ## 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 var.gamma <- rgamma(n=100, shape=2.5, scale=1) ## distribución gamma ## 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.333, 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.gamma) class(var.binomial3) class(varianza.normal) ## 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(3,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.gamma, density=5, freq=FALSE, main="distribución GAMMA", 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", breaks=10) 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 repetir el proceso con un N deseado mi_N <- 50 for (i in 1:10) { var.normal.mia <- rnorm(n=mi_N, mean=2.5, sd=1) par(mfcol=c(1,1)) par(mfcol=c(1,2)) hist(var.normal.mia, density=5, freq=FALSE, main="distribución simulada", ylim=c(0,0.8)) curve(dnorm(x, mean=mean(var.normal.mia), sd=sd(var.normal.mia)), col="red", lwd=2, add=TRUE, yaxt="n") qqnorm(var.normal.mia) qqline(var.normal.mia, 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 "fabricar" distribuciones de Poisson sindecimales <- as.integer(var.uniforme) ## quitando los decimales class(sindecimales) class(var.uniforme) sindecimales2 <- as.integer(round(var.uniforme, 0)) ## redondeando previamente el valor de la variable ## 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) ## hacemos nuestra primera matriz de datos (DATAFRAME) ## deben tener la misma longitud TODOS los vectores datos <- data.frame(var.uniforme, var.normal, var.poisson, var.gamma, var.negbin, var.binomial2, var.binomial3, mi.factor2, mi.factor3) ## 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$orden.azar <- runif(n=dim(datos)[1], min=0, max=1) datos.orden_azar <- datos[order(datos$orden.azar),] ## cambio los nombres de las variables names(datos.orden_azar) <- c("uniforme", "normal", "poisson", "gamma", "binomial_negativa", "binomial_01", "trinomial_0123", "factor_2", "factor_3", "azar") names(datos.orden_azar) ## renombramos el dataframe datos.originales <- datos datos <- datos.orden_azar ## elimino una variable del dataframe names(datos) datos$azar <- NULL names(datos) ## para seleccionar-extraer una parte de mis casos-observaciones y variables datos2 <- datos[c(1:15, 22:30),] datos3 <- datos[,c(1:6)] ## 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) ################################################################################ ## NUESTRO PRIMER MODELO FACTORIAL GENERAL ## ## una "NO-NORMAL" que no nos importa que no lo sea ... ¡¡Y LO HAREMOS BIEN!! ## ################################################################################ rm(list=ls()) cat("\014") dev.off() ## fijamos el proceso de aleatorización para que a todos nos salga igual set.seed(666) ## ## generamos valores para la variable respuesta según una normal v1 <- rnorm(n=10000, mean=10, sd=1.2) v2 <- rnorm(n=10000, mean=15, sd=1.2) ## asignamos códigos para los niveles (0 y 1) de un facor s1 <- rep(0, times=10000) s2 <- rep(1, times=10000) ## combinamos los valores creando las variables respuesta y factor respuesta <- c(v1, v2) factor <- c(s1, s2) class(factor) factor <- as.factor(factor) ## vemos la "no normalidad" de la variable respuesta hist(respuesta, breaks=25) qqnorm(respuesta); qqline(respuesta, col="red") ## valoramos cómo de distintos son los niveles 0 y 1 del factor ## con muestreo aleatorio de diferencias diferencia_1vs0 <- v2-v1 hist(diferencia_1vs0, breaks=25) ## la hipótesis nula es que la media de las diferencias de factor1-factor0 = 0 abline(v=0, col="blue", lwd=3) qqnorm(diferencia_1vs0); qqline(diferencia_1vs0, col="red") ## mediante percentiles vemos los intervalos de confianza de las diferencias factor1-factor0 quantile(diferencia_1vs0, probs=c(0.025, 0.975)) ## a alpha=0.05 con dos colas quantile(diferencia_1vs0, probs=c(0.005, 0.995)) ## a alpha=0.01 con dos colas quantile(diferencia_1vs0, probs=c(0.0005, 0.9995)) ## a alpha=0.001 con dos colas ## y lo visualizamos hist(diferencia_1vs0, breaks=25, main="distribución de diferencias factor1-factor0", sub="intervalo al 95% (alpha=0.05)") abline(v=quantile(diferencia_1vs0, probs=0.025), col="red", lwd=3) abline(v=quantile(diferencia_1vs0, probs=0.975), col="red", lwd=3) abline(v=0, col="blue", lwd=3) ## mediante un ANOVA "masivo" ## ya explicaremos la siguiente línea de código options(contrasts=c(factor="contr.sum", ordered="contr.poly")) ## modeloANOVA <- lm(respuesta ~ factor) summary(modeloANOVA) ## obtenemos los residuos del modeloANOVA y los estudiamos gráficamente hist(residuals(modeloANOVA), breaks=25) qqnorm(residuals(modeloANOVA)) qqline(residuals(modeloANOVA), col="red", lwd=1) ## otra forma equivalente a los normal probability plots de muchos paquetes estadísticos qqnorm(residuals(modeloANOVA), datax=TRUE) qqline(residuals(modeloANOVA), datax=TRUE, col="red", lwd=1) ## repetimos el ANOVA pero con una muestra de nuestros datos extraída al azar ## seleccionando solamente 25 datos por nivel del factor respuesta50 <- respuesta[c(101:125, 10101:10125)] factor50 <- factor[c(101:125, 10101:10125)] ## modeloANOVA50 <- lm(respuesta50 ~ factor50) summary(modeloANOVA50) ## obtenemos los residuos del modeloANOVA y los estudiamos gráficamente hist(respuesta50, breaks=10) hist(residuals(modeloANOVA50), breaks=10) qqnorm(respuesta50, datax=TRUE) qqline(respuesta50, datax=TRUE, col="red", lwd=1) qqnorm(residuals(modeloANOVA50), datax=TRUE) qqline(residuals(modeloANOVA50), datax=TRUE, col="red", lwd=1) ## revertimos el haber fijado la extracción aleatoria de datos set.seed(NULL)