p = seq(0.01, 0.99, by=0.01)
fsqrt(p)
max(abs(fsqrt(fsqrt(p), inverse=TRUE) - p)) # Should be 0
p = c(seq(-0.02, 0.02, by=0.01), seq(0.97, 1.02, by=0.01))
fsqrt(p)  # Has NAs
p = seq(0.01, 0.99, by=0.01)
par(mfrow=c(2,2))
y = seq(-4, 4, length=100)
for(d in 0:1) {
    matplot(p, cbind(logit(p, deriv=d), fsqrt(p, deriv=d)),
            type="n", col="purple", ylab="transformation",
            lwd=2, las=1,
            main=if(d==0) "Some probability link functions"
            else "First derivative")
    lines(p, logit(p, deriv=d), col="limegreen", lwd=2)
    lines(p, probit(p, deriv=d), col="purple", lwd=2)
    lines(p, cloglog(p, deriv=d), col="chocolate", lwd=2)
    lines(p, fsqrt(p, deriv=d), col="tan", lwd=2)
    if(d==0) {
        abline(v=0.5, h=0, lty="dashed")
        legend(0, 4.5, c("logit", "probit", "cloglog", "fsqrt"),
               col=c("limegreen","purple","chocolate", "tan"), lwd=2)
    } else
        abline(v=0.5, lty="dashed")
}
for(d in 0) {
    matplot(y, cbind(logit(y, deriv=d, inverse=TRUE),
                     fsqrt(y, deriv=d, inverse=TRUE)),
            type="n", col="purple", xlab="transformation", ylab="p",
            lwd=2, las=1,
            main=if(d==0) "Some inverse probability link functions"
            else "First derivative")
    lines(y, logit(y, deriv=d, inverse=TRUE), col="limegreen", lwd=2)
    lines(y, probit(y, deriv=d, inverse=TRUE), col="purple", lwd=2)
    lines(y, cloglog(y, deriv=d, inverse=TRUE), col="chocolate", lwd=2)
    lines(y, fsqrt(y, deriv=d, inverse=TRUE), col="tan", lwd=2)
    if(d==0) {
        abline(h=0.5, v=0, lty="dashed")
        legend(-4, 1, c("logit", "probit", "cloglog", "fsqrt"),
               col=c("limegreen","purple","chocolate", "tan"), lwd=2)
    }
}
# This is lucky to converge
earg = list(min=0, max=1, mux=5)
data(hunua)
fit.h = vglm(agaaus ~ bs(altitude),
             fam= binomialff(link="fsqrt", earg=earg),
             data=hunua, trace=TRUE, crit="d")
plotvgam(fit.h, se=TRUE, lcol="red", scol="red",
     main="Red is Hunua, Blue is Waitakere")
predict(fit.h, hunua, type="response")[1:3]
# The following fails.
data(pneumo)
pneumo = transform(pneumo, let=log(exposure.time))
earg = list(min=0, max=1, mux=10)
fit = vglm(cbind(normal, mild, severe) ~ let,
           cumulative(link="fsqrt", earg=earg, par=TRUE, rev=TRUE),
           data = pneumo, trace=TRUE, maxit=200)Run the code above in your browser using DataLab