Comparative graph of exponential models

Question

Comparative graph of exponential models

Navigation

#1 by (1 votes)

4

I have a 3x3 factorial experimental design and I want to perform a non-linear regression to adjust exponentially. I can easily adjust the curves individually in the following way:

dataUrea0 <- subset(Suelo, Nurea == "0")
dataUrea90 <- subset(Suelo, Nurea == "90")
dataUrea180 <- subset(Suelo, Nurea == "180")

plotPoints(P ~ NCP, data = dataUrea0)
nlsUrea0 <- nls(P~a*exp(b*NCP), dataUrea0, start=list(a= 5, b=0.04))
summary(nlsUrea0)
pcrGOF(nlsUrea0, PRESS = FALSE)
overview(nlsUrea0)
plotfit(nlsUrea0, smooth = TRUE, xlab="Norg (Kg/ha)", ylab="P (ppm)", col.fit = "blue", lwd = 3)

plotPoints(P ~ NCP, data = dataUrea90)
nlsUrea90 <- nls(P~a*exp(b*NCP), dataUrea90, start=list(a= 5, b=0.04))
summary(nlsUrea90)
pcrGOF(nlsUrea90, PRESS = FALSE)
overview(nlsUrea90)
plotfit(nlsUrea90, smooth = TRUE, xlab="Norg (Kg/ha)", ylab="P (ppm)", col.fit = "blue", lwd = 3)

plotPoints(P ~ NCP, data = dataUrea180)
nlsUrea180 <- nls(P~a*exp(b*NCP), dataUrea180, start=list(a= 5, b=0.04))
summary(nlsUrea180)
pcrGOF(nlsUrea180, PRESS = FALSE)
overview(nlsUrea180)
plotfit(nlsUrea180, smooth = TRUE, xlab="Norg (Kg/ha)", ylab="P (ppm)", col.fit = "blue", lwd = 3)

However, when I want to make a comparative graph I get an error. I write in the r:

S <- within(S, {
  Nmin <- as.factor(NM)
})
S <- within(S, {
  Norg <- as.factor(NCP)
})
S$Tratamiento <- as.factor(S$Tratamiento)

S%>%
  split(.$NM) %>%
  map( ~nls(P~a*exp(b*Norg)),
              start = list(a = 5),
                           (b = 0.04),
              trace = TRUE, 
              algorithm = "port", 
              data=.) %>% 
  map_df(~augment(.), .id="NM") %>% 
  as.tibble() %>%                      
  ggplot(aes(x=Norg, y=P, color=NM)) +
  geom_line(aes(y=.fitted))+
  labs(title="",
       x="x",  
       y="y")+
  theme_minimal()

And I get the following error:

Error: '.x' is not a vector (language)

To do all this I use the libraries:

library(mosaic)
library(ggplot2)
library(nlstools)
library(minpack.lm)
library(qpcR)
library(broom)
library(purrr)
library(tidyverse)

Here is a sample of the data: link Thanks!

r

asked by germanfernandez 01.02.2018 в 04:12

source

1 answer

Upload an image in HTML relative paths Spring Boot with server.contextPath

score 1 · Accepted Answer

Germán,

Your code has two problems, one syntax and one data type problem.

The syntax problem is that there are several parentheses wrongly located in the line of map , that's why the function does not find the data and returns the error .x is not a vector . In these cases I suggest you try the code that you are going to vectorize with map() separately, that facilitates getting informative errors. When using map() the errors that come back are not very useful to find the problem. Also, since the syntax is simpler, it is easier to count the parentheses and verify the location of the commas, which is what is failing in your code.

If the syntax error is solved, a second error will appear: data type: Norg should not be a factor, but numerical. If you look at your code you will see that you use them for an arithmetic operation in b*Norg . Possibly nls() coercion that factor to numeric, giving integers greater than 1. It is possible that the model fails for that reason.

Solution

S%>%                                  # S son los datos leídos del xlsx en DropBox
  mutate(Norg=NCP) %>%                #Reemplaza a S <- within(S, {Norg <- as.factor(NCP)}) y además no lo hace factor.
  split(.$NM) %>%
  map( ~nls(P~a*exp(b*Norg),
       start = list(a = 5, b = 0.04),  # Acá había un error, (b = 0.04) quedaba fuera de la lista start
       trace = TRUE, 
       algorithm = "port", 
       data=.)) %>%                   #Acá faltaba un paréntesis.
  map_df(~augment(.), .id="NM") %>%   #De acá salen los warnigns
  as.tibble() %>%       
  mutate(NM=factor(NM, c("0", "90", "180"))) %>% #Paso NM a factor y le doy un orden específico, que se reflejará en las etiquetas de leyenda. 
  ggplot(aes(x=Norg, y=P, color=NM)) +
    geom_line(aes(y=.fitted))+
    labs(title = "",
             x = "x",  
             y = "y") +
    theme_minimal()

The graphic you're looking for returns to me, with a line for NM level, ordered from 0 to 180. This code works, although it produces warning that are not worrisome. I think that because of the use that augment() makes of an obsolete function. Nothing serious, for now.