0

Tengo varios datos los cuales algunos los modelo como una regresión no lineal (nlsLM) sigmoidal y otros como una regresión lineal (lm) cuadrática y necesito representarlas en el mismo gráfico. Para esto, necesito poner los datos en porcentaje absoluto (que cada columna sea expresada en porcentaje siendo el 100% el valor máximo de cada columna), por tener escalas distintas. Sé que esto se puede realizar fácilmente en excel, pero también me gustaría aprender a hacerlo en R. Para los datos de M, N, P y K yo los ajusto de la siguiente manera:

    plotPoints(N ~ Dia, data = data)
NT1 <- nlsLM(N ~ a + ((b - a)/(1 + exp(-c * (Dia - d)))),
              data = data,
              start = list(a = min(data$N),
                           b = max(data$N),
                           c = 1, d = median(data$Dia)),
              trace = TRUE, algorithm = "port")
pcrGOF(NT1, PRESS = FALSE)
overview(NT1)
plotfit(NT1, smooth = TRUE, xlab="Tiempo (dds)", ylab="KgN/ha", col.fit = "blue", lwd = 3)

Y para el parametro I lo realizo de la siguiente manera:

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

El dataframe es lo dejo aquí: https://www.dropbox.com/s/g9ic9e0vcau9cvu/data.xlsx?dl=0 Se observa también la dificultad que varios de los datos son N/A, para evitar esto use la función predict para predecir los valores en blanco, estos están en las columnas variable_predict_vals.

Muchas gracias!

5
  • Hola en que paquetes están la funciones que usas? Commented el 11 abr. 2018 a las 18:43
  • Tampoco se puede reproducir el primer gráfico, plotPoints, el error Error in .plotMaps(polys = data, xlim = xlim, ylim = ylim, projection = projection, : Invalid PolySet. Hay algún proceso que hagas sobre los datos? Commented el 11 abr. 2018 a las 21:00
  • Tampoco se puede ajustar el modelo. En este caso el problema son los NA, posiblemente para el gráfico también. En la respuesta imputé 0 a los missing para probar el código y así se soluciona la parte técnica. Los resultados no son nada confiables.
    – mpaladino
    Commented el 11 abr. 2018 a las 23:26
  • Hola, uso estas librerías: library(mosaic) library(ggplot2) library(nlstools) library(minpack.lm) library(qpcR) library(broom) library(purrr) library(tidyverse) library(scales) Puede que el gráfico de error porque cuando los gráfico individualmente borro las columnas que contienen los N/A, estos están por problemas en el laboratorio. Este problema lo solucione recién con la función predict. Commented el 12 abr. 2018 a las 15:37
  • Va respuesta nueva al problema técnico del sacar el gráfico. Insisto en que el problema de manejo de missing es muy complicado. Si tenés seguridad de que tus missing son completamente aleatorios podrías manejarlos "trivialmente": eliminarlos, imputarlos por modelo. Pero deberían revisar a fondo los protocolos y procesos de laboratorio, si esto fuera para publicación. Si los missing no son completamente aleatorios podrían comprometer los resultados de investigación. dev.off(profe_de_metodologia)
    – mpaladino
    Commented el 13 abr. 2018 a las 3:23

1 respuesta 1

0

Con los datos "completos" el gráfico se puede hacer.

Esta respuesta en SO https://stackoverflow.com/questions/45971661/how-to-use-method-nlslm-in-packages-minpack-lm-in-geom-smooth (que te había comentado hace un tiempo) tiene el código básico para resolver el problema. Acá está adaptado para varios modelos.

library(readxl)        # Por qué no dput()?
library(tidyverse)
library(scales)        # Para rescale(). 
library(minpack.lm)    # Para ajustar nlsLM.

datos <- read_excel("data.xlsx")

Los datos son los del DropBox. Llamar data a los datos no es una buena práctica. Si bien no es una palabra reservada de R forma parte de unas cuantas. No es buena práctica llamar a los datos data.

Graficar varios modelos.

Si lo que estás buscando es graficar las curvas de ajuste de cada modelo (1 por cada variable) se puede hacer usando color = dentro de ggplot(aes()). En este caso se grafican 5 curvas por cada tipo de modelo y son dos tipos de modelos, el nlmLM u un OLS exponencial.

library(readxl)        # Por qué no dput()?
library(tidyverse)
library(scales)        # Para scales. 
library(minpack.lm)    # Para ajustar nlsLM.

datos <- read_excel("./gustavo/data.xlsx")

datos %>% 
  select(Bloque, Dia, ends_with("vals")) %>%       #Elimino las variables que sobran. 
  mutate_at(vars(ends_with("vals")), rescale) %>%  #Reescalo todas las numericas a un valor entre 0 y 1. 
  gather(variable, valor, -Bloque, -Dia) ->        #Paso los datos a formato largo, el subset lo hará color= en ggplot()
  datos_largos  

ggplot(datos_largos, aes(x = Dia, y= valor, color = variable), ) +
  geom_point() +
  geom_smooth(method = "nlsLM", 
              formula = y ~ a + ((b - a)/(1 + exp(-c * (x - d)))),  #En la fórmula uso los nombres del variable definidos en ggplot(aes())
              # Comentario a la línea siguiente:
              # Acá hay un posible problema si no se conocen previamente los parámetros de inicio del modelo.
              # En este caso funciona porque los datos están estandarizados 0-1, entonces sé de antemano máximo y mínimo. Pero...
              # No encuentro una manera de que ggplot calcule estos parámetros en runtime para cada subset definido por color =. 
              # En la otra respuesta se puede hacer pidiendo las variables en el entorno global. min(datos$N), datos está en Global_environment
              # En este caso no se puede porque los datos están en formato largo y el subset lo hace ad hoc ggplot(aes(color = ))
              # No sé cómo hacer que stat_smooth() vaya al parent.frame() de ggplot(), donde sí están y y x y color controla el subset. 
              # Va a buscarlos a directamente Global_environment. 
              method.args = list(start = list(a = 0,      
                                              b = 1,
                                              c = 1, 
                                              d = 58.87)
                                 trace = TRUE, 
                                 algorithm = "port"),
              se = FALSE) +                               # Importante: si TRUE, falla.
  geom_smooth(formula = y ~ x + I(x^2)-1, 
              method = "lm",
              se=FALSE, 
              linetype= 5) +                              # Guiones para el OLS exponencial. No sale en etiquetas. 
  theme_minimal()

Resultado:

introducir la descripción de la imagen aquí

Funciona siempre que conozcas de antemano los parámetros de inicio de nlsLM y que sean los mismos para todos los modelos. No encuentro una manera de forzar a stat_smooth() -que es el que hace el trabajo pesado de geom_smooth()- a que busque las variables en el entorno de la función que lo llama. Las busca en el entorno global.

4
  • Estimado, el problema de los N/A lo resolví usando la función predict, pero me encontré ante los siguientes problemas: - Cuando trato de gráficar M, N, P, K e I juntos, copio lo que escribiste con N y lo pego a continuación remplazando por M, P, K e I, y juntos me saltan varios errores, -En el eje de las "x" necesito que aparezcan los datos de los días en que fueron tomadas las muestras y no el porcentaje absoluto (eso es solo para el eje de las "y", perdón si no fui claro), -Al tener errores no pude probar, pero necesitaría una leyenda en el gráfico. Muchas gracias! Commented el 12 abr. 2018 a las 15:23
  • Hola, no entiendo muy bien. Lo que buscas es una curva de ajuste del nlsLM para cada variable (N, M, P, K) y, aparte, ajustar para I un modelo lineal OLS? O ajustar ambos modelos a todas las variables, de N hasta I? En este último caso tendrías 10 curvas, 2 modelos en 5 variables.
    – mpaladino
    Commented el 12 abr. 2018 a las 19:00
  • Podrías cargar los datos con el problema de los missing resuelto? Con los 0 que agregué los modelos no convergen y dan error, interrumpiendo la ejecución. Entonces no puedo saber si lo que falla es el código o los datos.
    – mpaladino
    Commented el 12 abr. 2018 a las 21:24
  • Ya edite el documento de excel agregándole las columnas con los valores predecidos y edite la pregunta. Disculpa si mi omisión te causo muchas dificultades. Gracias Commented el 12 abr. 2018 a las 23:56

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.