Ejemplo
1 enunciado
El fichero adjunto contiene información sobre animales bovinos de la raza Asturiana de los Valles. Algunas variables son:
- IA
- inseminación artificial; si se trata de un semental.
- manejo
- alimentación (1 = sólo leche; 2X = leche + piensoX)
- peso
- peso al destete.
- calificacion
- morfológica.
- conformacion
- de la canal (SEUROP = S/EUROP) S=superior.
- numero
- de parto (lactación).
- facilidad
- dificultad al parto (1 = pare sola ; 2 = ligera ayuda ; 3 = fuerte tracción ; 4 = cesárea).
- tipo
- normal, aculonado, culón (a ojo)
- genotipo
- normal, heterocigoto, culón (laboratorio)
Tareas:
- Calcular
- edadd
- edad al destete
- edads
- edad al sacrificio
- Determinar automáticamente qué variables son cualitativas y cuáles cuantitativas.
- Buscar un criterio para eliminar (sustituir por NA) los valores erróneos.
- Automáticamente, elimina los valores negativos de las variables cuantitativas y el 1% de los datos extremos (0,5% superior y 0,5% inferior).
- Analiza la relación entre
- tipo y genotipo.
- calificacion y: alzada, longcruz, anchomuslos, curvanalga, peritorax, longgrupa.
- peso y edadd.
- Busca cómo predecir a partir del resto de variables:
- genotipo
- facilidad
2 posible resolución (parcial)
load("vacas.rda") # el dataframe se llama "v" V <- v # copiarlo para conservar los originales ## 1. Calcular ## - edadd :: edad al destete ## - edads :: edad al sacrificio fecha <- function (x) as.Date(x, "%d/%m/%Y") V$f_destete <- fecha(v$f_destete) V$f_nacim <- fecha(v$f_nacim) V$FSacrificio <- fecha(v$FSacrificio) V$edadd <- as.numeric(V$f_destete - V$f_nacim) # as.numeric para evitar rarezas de difftime V$edads <- as.numeric(V$FSacrificio - V$f_nacim) ## 2. Determinar automáticamente qué variables son cualitativas y cuáles ## cuantitativas. ## cualitativa = character o factor ## cuantitativa = numeric o Date es.cualitativa <- function (x) is.character(x) | is.factor(x) es.cuantitativa <- function (x) is.numeric(x) | inherits(x,"Date") cualitativas <- sapply (V, es.cualitativa) cuantitativas <- sapply (V, es.cuantitativa) ## comprobación stopifnot(all(xor(cualitativas, cuantitativas))) stopifnot(sum(cuantitativas)+sum(cualitativas) == ncol(V)) # lo mismo ## otra forma, para evitar la comprobación: es.cuantitativa <- Negate(es.cualitativa) ## 3. Buscar un criterio para eliminar (sustituir por NA) los valores ## erróneos. summary (V[cualitativas]) ## parece que no hay niveles erróneos en las cualitativas; ## comprobar nomMunicipio porque tiene muchos: table(V$nomMunicipio) # sobran muchos niveles V$nomMunicipio <- droplevels(V$nomMunicipio) ## en el tipo, AC debería estar entre C y N: V$tipo <- factor (v$tipo, levels=c("C ","AC","N ")) ## la conformación está mal ordenada: table(v$conformacion) levels(v$conformacion) ## de los niveles SEUROP indicados en el enunciado (ojo a los espacios) ## aparecen sólo SEUR; además, cada letra puede estar seguida de un signo menos; ## en el orden adecuado S sería el nivel superior y R el inferior: niveles.conf <- c("R ","U-","U ","E-","E ","S ") V$conformacion <- factor (v$conformacion, levels=niveles.conf) ## no se aprecian errores en las cualitativas; summary(V[cuantitativas]) ## los errores de las cuantitativas se tratan en el apartado siguiente ## 4. Automáticamente, elimina los valores negativos de las variables ## cuantitativas y el 1% de los datos extremos (0,5% superior y 0,5% inferior). eliminar.negativos <- function (x) {x[x<=0] <- NA; x} ## type=2 para evitar interpolar entre fechas: eliminar.extremos <- function (x) {u <- quantile(x,c(.005,.995),na.rm=TRUE,type=2); x[x<u[1] | x>u[2]] <- NA; x} eliminar.raros <- function (x) eliminar.extremos (if (is.numeric(x)) eliminar.negativos(x) else x) for (i in which(cuantitativas)) V[,i] <- eliminar.raros(V[,i]) summary (V[cuantitativas]) pdf("vcuan.pdf") for (i in which(cuantitativas)) if (length (unique (V[,i])) < 10) barplot (table (V[,i]), main=names(V)[i]) else hist(V[,i], breaks=10, main=names(V)[i]) dev.off() ## las distribuciones parecen razonables ## 5. Analiza la relación entre ## - tipo y genotipo. ## - calificacion y: alzada, long_cruz, ancho_muslos, curva_nalga, peri_torax, ## long_grupa. ## - peso y edadd. table(V$genotipo, V$tipo) ## los heterocigotos se clasifican a ojo casi todos como normales ## (el alelo normal es recesivo: mh/+ se comporta como +/+) modelo <- lm (calificacion ~ alzada + long_cruz + ancho_muslos + curva_nalga + peri_torax + long_grupa, V) summary(modelo) # explica solamente el 18% ## para simplificarlo: ##step(modelo) # da error por los NA regresores <- c("alzada", "long_cruz", "ancho_muslos", "curva_nalga", "peri_torax", "long_grupa") Vmodelo <- na.omit (V[c("calificacion",regresores)]) formula <- paste ("calificacion ~", paste (regresores, collapse="+")) modelo <- lm(formula, Vmodelo) summary(modelo) # R2 = 18% modelo.simple <- step(modelo) summary(modelo.simple) # long_grupa discutible summary(lm(calificacion~long_cruz+ancho_muslos+curva_nalga+peri_torax,V)) #R2=18% ## parece que la definición de calificación ha cambiado con el tiempo: pdf("vcalif.pdf") plot (calificacion ~ f_nacim, V) plot (calificacion ~ f_nacim, V, subset=f_nacim>as.Date("2000-01-01")) dev.off() ## peso al destete en función de la edad al destete: summary (m1 <- lm (peso ~ edadd, V)) # medio kilo diario; R2=31% summary (m2 <- lm (peso ~ edadd + I(edadd^2), V)) # R2=32% ## el término cuadrático es significativo pero aporta poco al R2 pdf("vpeso.pdf") plot (peso ~ edadd, V) abline (m1, col=2, lwd=3) x <- sort(V$edadd) y <- predict (m2, data.frame(edadd=x)) lines (x, y, col=3, lwd=3) # tiene sentido el cuadrático dev.off() ## 6. Busca cómo predecir a partir del resto de variables: ## - genotipo ## - facilidad round (sapply(V[cualitativas], function(var) chisq.test(V$genotipo, var)$p.value), 3) ## => descartar IA como predictor de genotipo round(sapply(V[cuantitativas], function(var) kruskal.test(var,V$genotipo)$p.value), 3) ## => descartar pesos y variables relacionadas con fecha de sacrificio pdf("vgenotipo.pdf") for (i in which(cuantitativas)) boxplot(as.formula(paste(names(V)[i],"~genotipo")), V, notch=TRUE, main=names(V)[i]) dev.off()