1
$\begingroup$

Topic: probit and logprobit regression.

Context: I have to implement a model of size at maturity in a wild population and I must choose a binomial linear model to estimate the parameters based on imperfectly observed reproductive events. There are theoretical reasons to expect that size-at-maturity to be more log-normally than normally distributed. I am using R's nls() to scout models which will ultimately be implemented as part of a larger JAGS model. To account for imperfect detection, I need the model asymptote < 1 which is why I am using nls rather than glm() or similar.

Question: Why does y~plnorm(x*coef + Int) produce very different results than y~pnorm(log(x)*coef + Int) (which is closer to y~pnorm(x*coef + Int)?

Code to simulate data, fit models, plot 'data' and model fits:

dummy <- function(){
  # Dummy data
  detection <- 0.75
  size      <- (1:35)/2
  quants    <- 1.5*size - 18
  pvec      <- detection*plnorm(quants)
  pvec      <- as.matrix(pvec)
  n         <- 200
  obs       <- t(apply(pvec, 1, FUN=function(x){ rbinom(n,1,x) }))

  success   <- apply(obs, 1, FUN=sum)
  aggr.obs  <- data.frame(size, n, success)
  for(i in 1:dim(aggr.obs)[1]){
    bt              <- binom.test(x=aggr.obs[i,"success"], n=aggr.obs[i,"n"])
    aggr.obs$est[i] <- bt$estimate
    aggr.obs$lcl[i] <- bt$conf.int[1]
    aggr.obs$ucl[i] <- bt$conf.int[2]
  }

  colnames(obs)  <- paste("X", 1:n, sep="")
  obs.dum        <- data.frame(size=size, obs)
  obs.dum        <- reshape(obs.dum, varying=list(colnames(obs)), direction="long")
  obs.dum        <- obs.dum[,c(1,3)]
  names(obs.dum) <- c("size", "obs")
  return(list(aggr.obs=aggr.obs, obs.dum=obs.dum, pvec=pvec))
}
x        <- dummy()
dat.long <- x$obs.dum
dat.agg  <- x$aggr.obs

# probit
nls.p <- nls(obs~pnorm(size*coef + Int)*detect, data=dat.long,
             start=list(coef=0.38, Int=-5, detect=0.8))

# lognormal probit 1
nls.logp <- nls(obs~pnorm(log(size)*coef + Int)*detect, data=dat.long,
                start=list(coef=16, Int=-41, detect=0.75))
# lognormal probit 2
nls.logp2 <- nls(obs~plnorm(size*coef + Int)*detect, data=dat.long,
                 start=list(coef=1.5, Int=-17, detect=0.75))

newdata    <- data.frame(size=seq(from=0.5, to=18, length=100))
pred.p     <- predict(nls.p,     newdata=newdata)
pred.logp  <- predict(nls.logp,  newdata=newdata)
pred.logp2 <- predict(nls.logp2, newdata=newdata)

nls.pred <- rbind(
  data.frame(size=newdata$size, fit=pred.p,     model="probit"),
  data.frame(size=newdata$size, fit=pred.logp,  model="logprobit"),
  data.frame(size=newdata$size, fit=pred.logp2, model="logprobit2")
)

ggplot(data=nls.pred, aes(x=size, y=fit, colour=model)) +
  geom_line() +
  geom_point(data=dat.agg, aes(x=size, y=est, colour=NULL)) +
  geom_errorbar(data=dat.agg, aes(x=size, y=est, ymin=lcl, ymax=ucl, 
                                    colour=NULL), width=0)

Model fits vs simulated data

$\endgroup$
1
$\begingroup$

I believe have figured out the answer. Explicitly including the default variables passed to the pnorm and plnorm function calls makes it clearer:

# lognormal probit 1
nls.logp <- nls(obs~pnorm(log(size)*coef + Int, mean = 0, sd = 1)*detect, data=dat.long, start=list(coef=16, Int=-41, detect=0.75))

# lognormal probit 2
nls.logp2 <- nls(obs~plnorm(size*coef + Int, lmean = 0, lsd = 1)*detect, data=dat.long, start=list(coef=1.5, Int=-17, detect=0.75))

Probit models work using a standard normal link function (mean = 0, sd = 1) and fitting parameters (in this case a coefficient and an intercept) that effectively translate the independent variables in quantiles of that standard normal distribution. However, using a logmean of 0 and and a logsd of 1 produces a highly skewed lognormal distribution which is not the shape as a standard normal distribution of log(indep. variable) transformed back to normal scale. If I change the formula to fit the log mean and log sd directly:

# lognormal probit 2
nls.logp2 <- nls(repro~plnorm(size, mu, sd)*detect, data = frepro,
             start = list(sd = 0.06, detect=0.79, mu=2.54))

I get the exact same fitted curve as I did for lognormal probit 1 using nls(obs~pnorm(log(size)*coef + Int, mean = 0, sd = 1)*detect

EDIT: However, lognormal probit 1 is not as good a fit (the data were simulated under lognormal probit 2). In order to fit a flexible lognormal probit that can fit data simulated under either model I found that an extra parameter is needed. The following works for this:

    #logprobit 3
    d.logprobit3 = nls(obs~plnorm(size+Int, logmean, logsd)*detect, 
       data = dat.long, detect=0.75), 
       start = list(Int=-10, logmean = 1, logsd = 0.5, detect=0.75))

EDIT2: This model is basically the three parameter lognormal CDF, so Int should have been called parameter theta, the location parameter. Because I also need this model to work in JAGS I have tried to find a way to fit it using pnorm() so that I can use JAGS' probit() link function (there is no plnorm equivalent link function that I know of). The barrier is that useful values of theta in log(size - theta) result in trying to log a negative. However the following hack works and allows me to fit the three parameter model using pnorm:

nls.logp4<-nls(obs~pnorm(log(sapply(size,FUN = function(x){max(x-theta,0)}))
                         *coef + Int)*detect, data = dat.long,
               start = list(theta = 12, coef = 1, Int = 1, detect=detect))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.