r nls gnm

Fitting Generalized Nonlinear Model in R

I want to fit the following Generalized Nonlinear Model: Probit(G)=K+1/Sigma*(Temp-T0)*Time. As naive model, I created Probits(G) by qnorm(G) and then fitted the Nonlinear Model. But I want to fit Nonlinear Model with logit link similar to glm function in R.

How can I fit such Generalized Nonlinear Model with logit link in R?

Data <-
  structure(list(Temp = c(23L, 23L, 23L, 23L, 23L, 23L, 23L, 23L,
  27L, 27L, 27L, 27L, 27L, 27L, 33L, 33L, 33L, 33L, 33L, 33L, 33L,
  35L, 35L, 35L, 35L, 35L), Time = c(144L, 168L, 192L, 216L, 240L,
  264L, 288L, 312L, 120L, 144L, 168L, 192L, 216L, 240L, 72L, 96L,
  120L, 144L, 168L, 192L, 216L, 96L, 120L, 144L, 168L, 192L), G = c(15,
  25.5, 27, 28, 28.5, 39.5, 41.5, 43, 13, 21.5, 29.5, 30.5, 32.5,
  35, 13.5, 28, 32.5, 33.5, 35, 39.5, 42, 6.5, 30, 39.5, 57, 58.5
  )), .Names = c("Temp", "Time", "G"), class = "data.frame", row.names = c(NA,
  -26L))

Data$GermRate <- 1/Data$Time
Data$Probits <- qnorm(p=Data$G/100) # Get Probits


fm1 <-
  nls(
      formula= Probits ~ K+1/Sigma*(Temp-T0)*Time
    , data=Data
    , start=list(K=1, Sigma=2, T0=2)
    #, algorithm= "port"
    )
fm1Summary <- summary(fm1)
fm1Coef <- summary(fm1)$coef

Solution

You can fit this type of model using the gnm package for generalized nonlinear models. It takes a bit of work, as gnm uses pre-defined functions of class "nonlin" to specify nonlinear terms in the model and the ones provided by the package are generally insufficient to specify an arbitrary nonlinear function. However it is possible to define a custom "nonlin" function to use with gnm.

In your model, k is a linear parameter, so we only need to worry about the second term. This can be specified via the following "nonlin"function

customNonlin <- function(Temp, Time){
    list(predictors = list(sigma = 1, t0 = 1),
         variables = list(substitute(Temp), substitute(Time)),
         term = function(predLabels, varLabels) {
             sprintf("1/%s * (%s - %s) * %s",
                     predLabels[1], varLabels[1], 
                     predLabels[2], varLabels[2])
         })
}
class(customNonlin) <- "nonlin"

In the returned list,

predictors specifies that sigma and t0 are predictors with a single intercept term (i.e. are individual parameters).
variables specifies that there are two variables, to be supplied by the user via the Temp and Time arguments.
term specifies a function to create a deparsed mathematical expression of the term, given labels for the predictors and the variables.

More detail on "nonlin" functions can be found in Section 3.5 of the gnm vignette.

Now we can try fitting your model as follows

mod1 <- gnm(cbind(G, 100 - G) ~ customNonlin(Temp, Time), family = binomial, 
            data  = Data, start = c(1, 2, 2))

Note that, as in glm, an intercept is added to the formula by default, which here will represent k. Although the starting values are far from the solution, the gnm convergence criteria are met at this point, so the algorithm does not perform any iterations. In this case a better starting estimate of sigma is required for gnm to converge to something more sensible

mod2 <- gnm(cbind(G, 100 - G) ~ customNonlin(Temp, Time), family = binomial, 
            data  = Data, start = c(1, 2000, 2))

mod2

Call:
gnm(formula = cbind(G, 100 - G) ~ customNonlin(Temp, Time), family = binomial, 
    data = Data, start = c(1, 2000, 2))

Coefficients:
(Intercept)        sigma           t0  
     -2.589     1915.602        8.815  

Deviance:            53.53157 
Pearson chi-squared: 49.91347 
Residual df:         23

Actually it is possible to specify this model using the Mult function provided by gnm, as long as you don't mind re-parameterising the model:

mod3 <- gnm(cbind(G, 100 - G) ~ Mult(1, 1 + offset(Temp), offset(Time)), 
            family = binomial, data  = Data,
            start = c(1, 1/2000, -2))

mod3

Call:

gnm(formula = cbind(G, 100 - G) ~ Mult(1, offset(Temp) + 1, offset(Time)), 
    family = binomial, data = Data, start = c(1, 1/2000, -2))

Coefficients:
                             (Intercept)  
                               -2.588874  
Mult(., 1 + offset(Temp), offset(Time)).  
                                0.000522  
Mult(1, . + offset(Temp), offset(Time)).  
                               -8.815152  

Deviance:            53.53157 
Pearson chi-squared: 49.91347 
Residual df:         23

EDIT

The function of the parameters is specified in the term component of the list returned by customNonlin, which you can see via

customNonlin(Temp, Time)$term(c("sigma", "t0"), c("Temp", "Time"))
"1/sigma * (Temp - t0) * Time"

So if you just want to change the functional form, you will need to modify the term function. If you want to add/remove parameters, you would also need to modify the list in the predictors component. Similarly if the new term requires you to add/remove variables, you would modify the variables component.