# Difference between probit of log(x) and probit using inverse lognormal cdf of x

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)


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))