Ejemplo

Índice

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:

  1. Calcular
    edadd
    edad al destete
    edads
    edad al sacrificio
  2. Determinar automáticamente qué variables son cualitativas y cuáles cuantitativas.
  3. Buscar un criterio para eliminar (sustituir por NA) los valores erróneos.
  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).
  5. Analiza la relación entre
    • tipo y genotipo.
    • calificacion y: alzada, longcruz, anchomuslos, curvanalga, peritorax, longgrupa.
    • peso y edadd.
  6. 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()

Autor: Carlos Enrique Carleos Artime

Created: 2023-02-06 lun 17:32

Validate