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), lwd = (mylwd <- 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", las = 1,
          main = if (d == 0) "Some probability link functions"
          else "First derivative")
  lines(p, logit(p, deriv = d), col = "limegreen")
  lines(p, probit(p, deriv = d), col = "purple")
  lines(p, cloglog(p, deriv = d), col = "chocolate")
  lines(p, fsqrt(p, deriv = d), col = "tan")
  if (d == 0) {
    abline(v = 0.5, h = 0, lty = "dashed")
    legend(0, 4.5, c("logit", "probit", "cloglog", "fsqrt"), lwd = 2,
           col = c("limegreen","purple","chocolate", "tan"))
  } 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")
  lines(y, probit(y, deriv = d, inverse = TRUE), col = "purple")
  lines(y, cloglog(y, deriv = d, inverse = TRUE), col = "chocolate")
  lines(y, fsqrt(y, deriv = d, inverse = TRUE), col = "tan")
  if (d == 0) {
    abline(h = 0.5, v = 0, lty = "dashed")
    legend(-4, 1, c("logit", "probit", "cloglog", "fsqrt"), lwd = 2,
           col = c("limegreen","purple","chocolate", "tan"))
  }
}
par(lwd = 1)
# This is lucky to converge
fit.h <- vglm(agaaus ~ bs(altitude), binomialff(link = fsqrt(mux = 5)),
             data = hunua, trace = TRUE)
plotvgam(fit.h, se = TRUE, lcol = "orange", scol = "orange",
         main = "Orange is Hunua, Blue is Waitakere")
head(predict(fit.h, hunua, type = "response"))
# The following fails.
pneumo <- transform(pneumo, let = log(exposure.time))
fit <- vglm(cbind(normal, mild, severe) ~ let,
            cumulative(link = fsqrt(mux = 10), par = TRUE, rev = TRUE),
            data = pneumo, trace = TRUE, maxit = 200)Run the code above in your browser using DataLab