3

Tengo varios datos de distintos tratamientos a los cuales les hice una regresión lineal(lm) y en otros casos le hice una regresión no lineal (nls) según lo que indica la bibliografía. Quiero comparar gráficamente los nls entre si al igual que los lm. Dejo una muestra de datos y un script con gráficos.

library(mosaic)
library(ggplot2)
library(nlstools)
library(minpack.lm)
library(qpcR)

dataT1 <- subset(data, Trat == "T1")
dataT3 <- subset(data, Trat == "T3")
dataT7 <- subset(data, Trat == "T7")
plotPoints(MS ~ Dia, data = dataT1)

MST1 <- nlsLM(MS ~ a + ((b - a)/(1 + exp(-c * (Dia - d)))),
            data = dataT1,
            start = list(a = min(dataT1$MS),
                         b = max(dataT1$MS),
                         c = 1, d = median(dataT1$Dia)),
            trace = TRUE, algorithm = "port")

pcrGOF(MST1, PRESS = FALSE)
overview(MST1)
plotfit(MST1, smooth = TRUE, xlab="Tiempo (dds)", ylab="KgMS/ha", col.fit = "blue", lwd = 3)

Idem con T3 y T7 en el factor MS.

En el factor IAF:

IAFT1 <- lm(IAF ~ Dia + I(Dia^2)-1, data = dataT1)
summary(IAFT1)

dataT1 %>%
  ggplot(aes(x=Dia, y=IAF)) +
  geom_point() +
  geom_smooth(formula = y ~ x + I(x^2)-1, 
              method = "lm",
              se=FALSE) +
  labs(title="",
       x="Tiempo (dds)",  
       y="IAF (m2)", 
       caption="") + 
  theme_minimal()

Idem con T3 y T7 en el factor IAF.

Los datos en: https://www.dropbox.com/s/1budbw9awzjfc7t/Dudas.xlsx?dl=0

¿Cómo podría comparar todos los IAF y las MS? ¿Y el IAF y MS de un mísmo tratamiento? Gracias!

1 respuesta 1

4

Para lm

Para los lm es muy fácil. Simplemente en el código de ggplot que estás usando cambias la segunda línea para agregar un mapeado de color a la variable Trat.

#Librerías necesarios para los dos bloques de código. 

library(tidyverse)    #Para map(), map_df(), as_tibble(), ggplot()   
library(minpack.lm)   #Para nlsLM()
library(broom)        #Para augment()

dataT1 %>%   
  ggplot(aes(x=Dia, y=IAF, color=Trat)) +
  geom_point() +
  geom_smooth(formula = y ~ x + I(x^2)-1, 
              method = "lm",
              se=FALSE) +
      labs(x="Tiempo (dds)",  
      y="IAF (m2)", 
      color="Tratamiento") + 
  theme_minimal()

introducir la descripción de la imagen aquí

El resultado son tres curvas de ajuste con diferentes colores. Como lo mapeas directamente en la llamada a la función lo hereda el resto, por eso colorea los puntos de geom_point() y hace que geom_smooth() ajuste tres modelos, uno para cada nivel del factor Trat, y grafique las curvas con colores diferentes entre sí, pero los mismos que usa geom_point(). Ya sale todo etiquetado.

Para nlsLM

Para los nlsLM se puede con un método similar, te recomiendo visitar esta respuesta que trata exactamente ese problema. Si puedes leer inglés ahí está la solución.

Alternativamente podés generar los puntos que ajusta el modelo y plotearlos. Va el código, probado para nls. Hay bastante uso de iteración sobre listas con purrr::map() y familiares. Se me complicaría explicar en esta respuesta todo lo que hace, aunque traté de aclararlo en los muchos comentarios. Si conoces lapply() es muy similar, sólo cambia ligeramente la sintaxis y es más estricto con el tipo de datos de salida.

MST1 %>%              #Es el df que está en el enlace de DropBox
  split(.$Trat) %>%   #Armo una lista con 3 data.frame, uno por nivel de Trat
#map es similar a lapply, pasa una función a cada elemento de una lista y regresa una lista.
#Uso el . (punto) como referencia a los datos porque estoy usando tuberías %>% 
# El símbolo ~ indica a map que la siguiente expresión es una función, me ahorro declarar una función anónima como lo haría con lapply.
  map( ~nlsLM(MS ~ a + ((b - a)/(1 + exp(-c * (Dia - d)))),
              start = list(a = min(.$MS),
                           b = max(.$MS),
                           c = 1, 
                           d = median(.$Dia)),
              trace = TRUE, 
              algorithm = "port", 
              data=.)) %>% 
#augment() crea un data.frame al "aumentado" con los datos originales, la predicción del modelo y el error
  map_df(~augment(.), .id="Trat") %>%  #.id es un argumento para map_df(), indica el nombre de la columna clave, necesario para pasar de lista a data.frame  
  as_tibble() %>%                      
  ggplot(aes(x=Dia, y=MS, color=Trat)) +   #Acá empieza el gráfico propiamente dicho, lo anterior fue prepararle los datos. 
  geom_point() +
  geom_line(aes(y=.fitted))    #a x lo hereda, mapeo un nuevo y para la línea. 

introducir la descripción de la imagen aquí

Atención: el código para NLM no hace exactamente lo mismo que el código para lm. Solamente plotea los puntos en los que hay una estimación del valor de y, geom_line() se encarga de unirlos y crear la línea. Por eso no se ve como una curva y son visibles los saltos. Se puede corregir usando predict() con más data points de x, de ese modo se obtienen más puntos predichos de y en la gráfica se asemeja más a una curva. Aunque en sentido matemático estricto siguen siendo puntos discretos al tener suficientes puntos y usar la función de suavizado de geom_line() se ve como si fuera continua.

PD: si estás muy perdido con ggplot2 te paso un enlace a una guía básica que hice para mis alumnos.

7
  • Debo felicitarte por tu excelente guía. ¡Está muy completa!
    – aldo_tapia
    Commented el 28 nov. 2017 a las 10:54
  • Realmente muy completa la respuesta y me ayudo a despejar muchas dudas. El problema que estoy teniendo es en instalar el paquete "minpack.nls", ya que me sale el error: Warning in install.packages : package 'minpack.nls' is not available (for R version 3.4.0) Commented el 28 nov. 2017 a las 14:05
  • Porque no existe... es minpack.lm. Ya lo corregí.
    – mpaladino
    Commented el 28 nov. 2017 a las 15:01
  • Perdón por la demora en la respuesta. Ok, ahora me da el error de que no encuentra la función "map". Las librerías que tengo abiertas son: agricolae, mosaic, ggplot2, nlstools, minpack.lm, qpcR y broom. Puede ser que esta función pertenezca a otra librería que no sea alguna de estas? Commented el 30 nov. 2017 a las 16:56
  • Faltaba cargar purrr o mejor aún tidyverse que de paso también carga otras necesarias.
    – mpaladino
    Commented el 30 nov. 2017 a las 19:12

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.