Considera los datos¹ disponibles mediante SSH en carleos2.epv.uniovi.es:/home/manadine/dat/jacatón

## Tienes varias opciones para manejar los datos:
## - Trabajar remotamente usando por ejemplo
##   - el comando «ssh» 
##   - PuTTY
## - Descargar los ficheros en tu ordenador local mediante
##   - el comando «scp»
##   - el comando «sftp»
##   - Filezilla
##   - WinSCP

## Por ejemplo:

scp -r alumno@carleos2.epv.uniovi.es:/home/manadine/dat/jacatón ~/Documentos

Elabora un informe que incluya explicaciones someras sobre los siguientes pasos:

1. Describe los datos.

   • tamaños de los ficheros

     ## Si lo quieres hacer directamente en R:
     ped <- read.table ("~/Documentos/jacatón/pedigriHackathon.txt")
     met <- read.csv   ("~/Documentos/jacatón/MetadatosHackathon.csv")
     tig <- read.csv   ("~/Documentos/jacatón/ContigsTPM.csv")
     keg <- read.csv   ("~/Documentos/jacatón/keggTPM.csv", sep=";")
     keg <- read.csv2  ("~/Documentos/jacatón/keggTPM.csv", dec=".") # equivale a la anterior
     

     Si los ficheros son demasiado grandes, puedes explorarlos en Bash:

     cd ~/Documentos/jacatón
     head * | cut -c -100      # el kegg parece con sólo una línea
     file *                    # indica CR como salto de renglón (archivo de Mac)
     recode /CR.. keggTPM.csv  # cambia CR por LF; si no tienes "recode": cat keggTPM.csv | tr '\r' '\n' > keggTPM-LF.csv
     wc *                      # 3 ficheros con 341 renglones (uno de cabecera), 1 fichero con 2996
     awk '{print NF}' pedigriHackathon.txt  | uniq    # se ve a ojo
     awk -F ',' '{print NF}' ContigsTPM.csv | uniq
     awk -F ';' '{print NF}' keggTPM.csv    | uniq
     awk -F ',' '{print NF}' MetadatosHack* | uniq    # también a ojo
     ## Dividimos en bloques de 20000 variables para hacer más manejable el fichero de cóntigos:
     cut -d ',' -f -20000 ContigsTPM.csv > Contigs10000a.csv
     cut -d ',' -f 20001- ContigsTPM.csv > Contigs10000b.csv
     

   • tipos de variables y relaciones entre ellas
     • pedigrí
       • 3 variables con números enteros = identificadores de individuo
         • animal
         • padre
         • madre
       • el valor 0 indica dato desconocido
     • metadatos
       • x = identificador en el pedrigrí
       • reb = rebaño
       • RecordingYear = año
       • week = semana (junto con el año, darían la fecha de medición)
       • NLACTA = número de parto o lactación
       • DIASLE = número de días en la lactación
       • PGRASADC = porcentaje de grasa; variable respuesta
     • cóntigos
       • primera columna sin nombre = identificador del pedigrí
       • el resto de columnas son cuantificaciones de la presencia de cada cóntigo (tig000…)
     • keg
       • ID = identificador en el pedigrí
       • el resto de columnas son cuantificaciones de la presencia de cada función keg (K000…)

2. Consigue un único dataframe a partir de todas las tablas disponibles.

   mettigkeg <- cbind(met,tig,keg)
   todos <- merge (ped, mettigkeg, by.x = "V1", by.y = "x", all = TRUE)
   

3. Construye un modelo de predicción para PGRASADC.

   • explora variables individualmente

     r2 <- sapply (setdiff (names(todos), "PGRASADC"), #quito la variable dependiente 
                   function (variable) 
                     try (summary (lm (paste ("PGRASADC ~", variable), todos)) $ r.squared))
     mode(r2) <- "numeric" # permite conservar las etiquetas del vector
     

     • obtenemos un vector de \(R^2\)
     • usamos try para evitar que un error impida la creación del resultado
       • consecuencia: si hay algún error, el resultado deja de ser numérico y pasa a ser un vector de cadenas (character)
       • podríamos usar as.numeric para convertir a numérico, pero se perderían las etiquetas

       • podríamos recuperarlas así:

           R2 <- as.numeric(r2);  names(R2) <- names(r2)
           #+end_src R
         - para hacerlo de un solo paso, cambiamos directamene el modo: 
           #+begin_src R
           mode(r2) <- "numeric"
         

     • podemos intentar acelerar el proceso usando paralelismo
       • usaremos la función mclapply del paquete parallel
       • simplemente sustituyendo sapply por mclapply se trabaja con dos hilos de ejecucición (valor por omisión de la opción mc.cores)
       • detectCores() nos dice cuántos hilos en paralelo se podrían ejecutar

       • recomendación: dejar un hilo libre, para que el sistema responda sin problemas a la interacción con el usuario

         library (parallel)              # para mclapply
         nombres <- setdiff (names(todos), "PGRASADC")
         r2 <- as.numeric (mclapply (nombres,
                                     function (variable)
                                       try (summary (lm (paste ("PGRASADC ~", variable), todos)) $ r.squared),
                                     mc.cores = detectCores() - 1))
         names (r2) <- nombres           # mclapply no crea etiquetas a partir de su primer argumento
         

       • si necesitas medir el tiempo, puedes usar

         • system.time

           system.time ( {instrucción1 ; instrucción2 ; ...} )
           

         • proc.time

           comienzo <- proc.time()
           instrucción1 ; instrucción2 ; ...
           proc.time() - comienzo
           

     • exploramos los valores altos

       max (r2)                                 # hay un ajuste perfecto
       which.max (r2)                           # variable «Unclassified.»
       summary (todos$Unclassified.)            # factor pero debería ser numérica
       

       Vamos a intentar corregir la anomalía:

       niveles <- levels (todos$Unclassified.)
       niveles [is.na (as.numeric (niveles))    # qué valores no pueden convertirse en números
       ## vemos que se trata de valores con un punto al final;
       ## podemos pensar en varias opciones:
       ##  - eliminar el punto final
       ##  - quedarnos solamente con cuatro decimales
       ## vamos a desarrollar ambas
       
       ## eliminar punto final, versión 1:
       niveles.nuevos <- ifelse (is.na (as.numeric (niveles)), 
                                 substring (niveles, 1, nchar(niveles)-1),
                                 niveles)
       
       ## eliminar punto final, versión 2:
       niveles.nuevos <- gsub ("\\.$", "", niveles)
       
       ## mantener solo cuatro decimales, suponiendo que el punto decimal siempre está en la posición 7
       all (substring (niveles, 7, 7) == ".") # comprobación
       niveles.nuevos <- substring (niveles, 1, 11)
       
       ## mantener solo cuatro decimales, versión general 1:
       niveles.nuevos <- substring (niveles, 1, 4 + regexpr ("\\.", niveles))
       
       ## mantener solo cuatro decimales, versión general 2:
       niveles.nuevos <- gsub ("(.*\\..{4}).*", "\\1", niveles)
       
       ## en cualquier caso, hay que asignarlos
       levels (todos$Unclassified.) <- niveles.nuevos
       todos$inclasificables <- as.numeric (as.character (todos$Unclassified.))   # equivale a la siguiente
       todos$inclasificables <- as.numeric (niveles.nuevos) [todos$Unclassified.] # recomendado por ?factor
       r2[which.max(r2)]     <- summary (lm (PGRASADC ~ inclasificables, todos)) $ r.squared
       

   • construye un modelo con las mejores variables (explora cuántas puede asimilar tu ordenata)

     tail (sort(r2), 100)
     names (tail (sort(r2), 100))
     paste (names(tail(sort(r2),100)), collapse="+")
     paste ("PGRASADC~", paste(names(tail(sort(r2),100)),collapse="+"))
     lm (paste("PGRASADC~",paste(names(tail(sort(r2),100)),collapse="+")), todos)
     step (lm(paste("PGRASADC~",paste(names(tail(sort(r2),100)),collapse="+")),todos)) -> modelo
     

4. Describe la calidad del modelo.

   summary (modelo)
   

   Lo que nos interesa no es tanto la calidad del modelo en sí, sino la de la técnica para obtenerlo. Vamos a considerar solamente los registros con dato de PGRASADC.

   válidos <- todos [!is.na(todos$PGRASADC), ]
   

   • usa la misma técnica para calcular un modelo con el 80% de los datos

     entrena <- sample (nrow(válidos), 0.8*nrow(válidos)) # índices de los datos de entrenamiento
     modelo80 <- step (lm(paste("PGRASADC~",paste(names(tail(sort(r2),100)), collapse="+")),válidos[entrena,]))
     #                                                                                      ^^^^^^^^^^^^^^^^^
     summary (modelo80)  # calidad similar a «modelo»
     

   • estima el 20% restante

     estimas <- predict (modelo80, válidos[-entrena,])
     

   • compara con datos reales

     reales <- válidos[-entrena,"PGRASADC"]
     plot (reales ~ estimas)
     cor (estimas, reales) ^ 2
     summary (lm (reales ~ estimas))
     

   • repite eso varias veces y obtén valores medios

     correlas <- replicate (10,
                      {
                        entrena <- sample (nrow(válidos), 0.8*nrow(válidos)) # índices de los datos de entrenamiento
                        modelo80 <- step (lm(paste("PGRASADC~",paste(names(tail(sort(r2),100)), collapse="+")),válidos[entrena,]),
                                          trace = 0)
                        estimas <- predict (modelo80, válidos[-entrena,])
                        reales <- válidos[-entrena,"PGRASADC"]
                        cor (estimas, reales) ^ 2
                      })
     summary (correlas)