Nonlinear regression in R with dummy variables

Question

I need to perform three non-linear regression with the following formulas:

fire_vertical<-nls(R~0.23-exp(A*lnNt+C)+d*fire,start=list(A=-0.33,C=0.45,d=-0.001),trace=T)

fire_lateral<-nls(R~0.23-exp((A*lnNt)+(C+d*fire)),start=list(A=-0.33,C=0.45,d=-0.001),trace=T)

fire_nonlinear<-nls(R~0.23-exp((A+d*fire)*lnNt+C),start=list(A=-0.33,C=0.45,d=-0.001),trace=T)

I have the following data: R<-c(0.02,0.00,-0.06,0.11,0.06,0.00,-0.05,-0.06,0.02,-0.26,0.00,0.07,-0.07,0.23,0.06,-0.14,-0.04,0.09,-0.09,0.09,-0.02)

lnNt<-c(6.14,6.14,5.76,6.42,6.81,6.81,6.49,6.14,6.24,4.81,4.81,5.14,4.81,6.03,6.42,5.59,5.39,5.90,5.39,5.90,5.76)

fire<-c("before","before","before","before","before","before","before","afterone","afterone","afterone","afterone","afterone","afterone","afterone","aftertwo","aftertwo","aftertwo","aftertwo","aftertwo","aftertwo","aftertwo")

"lnNt" and "fire" are my independent variables. lnNt are continuous data, while "fire" is a categorical variable with three levels:"before", "afterone" and "aftertwo".

I need to transform my categorical variable ("fire") into a dummy variable to run the model. I'm not able to run the model with the dummy variable.

You could use package nlme, something like; `library(nlme); gnls(R ~ 0.79-exp(A*lnNt+C)+d, params = d ~ fire, data = DF, start = ...)`. Provide a reproducible example if you need further help. — Roland, Aug 05 '20 at 14:11
No, you provided a picture of your data (which is useless). Please read this: https://stackoverflow.com/a/5963610/1412059 — Roland, Aug 05 '20 at 14:47

Roland · Answer 1 · 2020-08-05T15:17:58.137

Use the gnls function from package nlme. It allows you to model parameters depending on a categorical variable (using the usual treatment contrasts by default).

DF <- data.frame(R, lnNt, fire, stringsAsFactors = TRUE)

fire_model<-nls(R~0.23-exp(A*lnNt+C)+d,
                start=list(A=-0.33,C=0.45,d=-0.001),trace=T, data = DF)
coef(fire_model)
#         A          C          d 
#-1.1367053  3.6708032 -0.1673437 

library(nlme)
fire_model_2 <- gnls(R~0.23-exp(A*lnNt+C)+d,params = list(A + C ~ 1, d ~ fire), 
                     start=c(coef(fire_model), rep(0, length(levels(DF$fire))- 1)), 
                     data = DF)

summary(fire_model_2)
# Generalized nonlinear least squares fit
# Model: R ~ 0.23 - exp(A * lnNt + C) + d 
# Data: DF 
#      AIC       BIC  logLik
# -33.3584 -27.09127 22.6792

# Coefficients:
#                    Value Std.Error    t-value p-value
# A              -1.247828  1.652336 -0.7551907  0.4611
# C               4.558928  7.565058  0.6026296  0.5552
# d.(Intercept)  -0.097251  0.128043 -0.7595203  0.4586
# d.fireaftertwo -0.062324  0.062210 -1.0018480  0.3313
# d.firebefore   -0.083946  0.064853 -1.2944170  0.2139
# 
# Correlation: 
#                A      C      d.(In) d.frft
# C              -0.999                     
# d.(Intercept)   0.864 -0.842              
# d.fireaftertwo  0.435 -0.454  0.050       
# d.firebefore   -0.148  0.117 -0.538  0.492
# 
# Standardized residuals:
#         Min          Q1         Med          Q3         Max 
# -1.66299973 -0.48743154  0.04367184  0.85175935  1.58042481

In this example we can clearly see that fire does not appear to have a significant impact.

fire_model <- gnls(R~0.23-exp(A*lnNt+C)+d,
                     start=coef(fire_model), data = DF)

anova(fire_model, fire_model_2)
#             Model df       AIC       BIC   logLik   Test  L.Ratio p-value
#fire_model       1  4 -35.05862 -30.88053 21.52931                        
#fire_model_2     2  6 -33.35840 -27.09127 22.67920 1 vs 2 2.299789  0.3167

Thank you so much for your answer. But actually, I'm doing model selection. This model that I initially put in the post is one of them (I have other variables too). For the fire, I want to test three models (I edited the post). Can the way you solve one of them apply to the others? Thank you. — Benjamin Guidon, Aug 05 '20 at 18:59
Model selection of non-linear models should always be based on domain knowledge, i.e. an understanding of the underlying processes. — Roland, Aug 06 '20 at 05:49
@Roland Can you please explain the syntax of `params = list(A + C ~ 1, d ~ fire)` ? What I am interpreting is: A and C from `start` parameter, and for `d ~ fire` estimate `d` based on the predictor `fire`. And when you say `fire` does not have a significant impact; this is based on the p-value of the `anova()` between the two models (one with the dummy variable, one without). Please correct my interpretation if not. — jack kelly, Dec 06 '22 at 05:35
`gnls` uses linear models for the dependence of parameters of the non-linear model on other variables. The syntax is explained in `help("gnls")`. The formula syntax for these is the same as for `lm`, only you can specify multiple parameters on the LHS. As a result, there are additional parameters in the non-linear model (as you can see in the summary output). The ANOVA output shows that including these additional parameters does not improve the model fit. — Roland, Dec 06 '22 at 06:23

Nonlinear regression in R with dummy variables

1 Answers1