Inverse probit

Whereas it is very simple to inverse the logit function [p=1/(1+exp(x)), x=log((1-p)/p)], I had some difficulties with the probit one. Here is the solution using pnorm().

> data <- data.frame(Doses=c(80, 120, 150, 150, 180, 200),

+                    Alive=c(10, 12, 8, 6, 2, 1),

+                    Dead=c(0, 1, 5, 6, 9, 15))

> ep <- glm(cbind(Alive, Dead) ~ Doses, data=data, family=binomial(link="probit"))

> ep$coefficients

(Intercept)       Doses 

 5.74329978 -0.03690572 

> ep$coefficients["(Intercept)"]+ep$coefficients["Doses"]*data$Doses

[1]  2.7908421  1.3146132  0.2074416  0.2074416 -0.8997301 -1.6378445

> pnorm(ep$coefficients["(Intercept)"]+ep$coefficients["Doses"]*data$Doses)

[1] 0.99737144 0.90568004 0.58216749 0.58216749 0.18413196 0.05072707

> predict(ep)

         1          2          3          4          5          6 

 2.7908421  1.3146132  0.2074416  0.2074416 -0.8997301 -1.6378445 

> predict(ep, type="response")

         1          2          3          4          5          6 

0.99737144 0.90568004 0.58216749 0.58216749 0.18413196 0.05072707