Search code examples
rnls

Comparing complete model versus a null model with nls

I'd like to comparing a complete non-linear model with a null model. Is this possible, considering approaches normally using with glm for example? In my case:

#Packages
library(minpack.lm)

# Data set - Diameter in function of Feature and Age
Feature<-sort(rep(c("A","B"),22))
Age<-c(60,72,88,96,27,
36,48,60,72,88,96,27,36,48,60,72,
88,96,27,36,48,60,27,27,36,48,60,
72,88,96,27,36,48,60,72,88,96,27,
36,48,60,72,88,96)
Diameter<-c(13.9,16.2,
19.1,19.3,4.7,6.7,9.6,11.2,13.1,15.3,
15.4,5.4,7,9.9,11.7,13.4,16.1,16.2,
5.9,8.3,12.3,14.5,2.3,5.2,6.2,8.6,9.3,
11.3,15.1,15.5,5,7,7.9,8.4,10.5,14,14,
4.1,4.9,6,6.7,7.7,8,8.2)
d<-dados <- data.frame(Feature,Age,Diameter)
str(d)

# Complet model
e1<- Diameter ~ a1 * Age^a2 
#Algoritm Levenberg-Marquardt
m1 <-  nlsLM(e1, data = d,
     start = list(a1 = 0.1, a2 = 10),
     control = nls.control(maxiter = 1000))

#Null model

e2<- Diameter ~ 1 
#Algoritm Levenberg-Marquardt
m0 <-  nlsLM(e1, data = d,
     control = nls.control(maxiter = 1000))
Warning message:
In nlsLM(e1, data = d, control = nls.control(maxiter = 1000)) :
  No starting values specified for some parameters.
Initializing ‘a1’, ‘a2’ to '1.'.
Consider specifying 'start' or using a selfStart model

Doesn't work and my final objective is:

anova(m1,m0)

Is this possible in non-linear universe?

Solution

  • Yes, there is an anova.nls defined in R

    methods("anova")
## [1] anova.glm*      anova.glmlist*  anova.lm*       anova.lmlist*  
## [5] anova.loess*    anova.mlm*      anova.nls*      anova.quantmod*
## see '?methods' for accessing help and source code

    however, the code in the question has some problems:

    • m1 should be using e2. As written both m0 and m1 are the same model.
    • e2 is incorrectly specified. It is specified in the question as if it were an lm model whereas it should be specified as an nls model.
    • I would be careful with using nlsLM. People often use it because it appears to converge when nls does not but that can occur because it prematurely converges giving an answer that is way off.

    Try this:

    o <- order(d$Age)

# Complet model
e1<- Diameter ~ Age^a2 
m1a <-  nls(e1, data = d[o, ], start = list(a2 = 10), alg = "plinear")

#Null model

ones <- rep(1, nrow(d))
e2 <- Diameter ~ a * ones
m0a <-  nls(e2, data = d[o, ], start = list(a = 1))

anova(m1a, m0a)
## Analysis of Variance Table
##
## Model 1: Diameter ~ Age^a2
## Model 2: Diameter ~ a * ones
##   Res.Df Res.Sum Sq Df  Sum Sq F value   Pr(>F)    
## 1     42     270.96                                
## 2     43     823.73 -1 -552.77   85.68 1.07e-11 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(Diameter ~ Age, d)
lines(fitted(m0a) ~ Age, d[o, ], col = "blue")
lines(fitted(m1a) ~ Age, d[o, ], col = "red")

    screenshot