[R] Differences between SPSS and R on probit analysis

[William Dunlap]

Did you not get a warning from glm, such as the following one?
> fm1 <- glm(affected/total ~ log(dose), family=binomial(link = probit), data=finney71[finney71$dose != 0, ])
Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
Do not ignore warnings.

The left hand side of the formula should a matrix containing the counts
of the afflicted and non-afflicted:
   cbind(affected, total-affected)
not the fraction of the total that were afflicted.  Then you would get

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -6.8940     1.0780  -6.395 1.60e-10 ***
log(dose)     0.9333     0.1344   6.944 3.82e-12 ***

Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Thu, Feb 23, 2017 at 12:26 PM, Biank M <bianca12_domi at hotmail.com> wrote:
> Hi,
>
> I'm working on the effects of alternative larvicides on Aedes aegypti. Right now, I am doing a binary mortality response with a single explanatory variable (dose) on 4 concentrations of one larvicide (+ control). Our university is fond of SPSS, and I have learned to conduct the basic probit model with it, including a natural logarithm transformation on my dosis data.
> Not so long ago, I've started working with R, and through a combination of the 'glm' and 'dose.p' functions, I get the same slope and intercept, as well as LD50 calculations. Nevertheless, the standard errors and Z-scores calculated through the Probit model in SPSS comes out completely different in R. Additionally, the 95% confidence intervals for the LD50 come out very differently between the two programs. I really don't have a clue on how I am getting the same slopes, intercepts and LD50's, but totally different SE, Z, and 95% CI. Can anybody help me so I can get the same results in R??
>
> I'll pass you the script and hypothetical data:
>
> dose <- c(6000, 4500, 3000, 1500, 0)
> total <- c(100, 100, 100, 100, 100)
> affected <- c(91, 82, 69, 49, 0)
>
> finney71 <- data.frame(dose, total, affected)
>
> fm1 <- glm(affected/total ~ log(dose),
> family=binomial(link = probit), data=finney71[finney71$dose != 0, ])
>
> xp1 <- dose.p(fm1, p=c(0.5,0.9))
> xp.ci <- xp1 + attr(xp1, "SE") %*% matrix(qnorm(1 - 0.05/2)*c(-1,1), nrow=1)
> EAUS.Aa <- exp(cbind(xp1, attr(xp1, "SE"), xp.ci[,1], xp.ci[,2]))
> dimnames(EAUS.Aa)[[2]] <- c("LD", "SE", "LCL","UCL")
>
> So, this is the regression results I get with R:
> summary(fm1)
>
> Deviance Residuals:
> 1 2 3 4
> 0.06655 -0.02814 -0.06268 0.03474
>
> Coefficients:
> Estimate Std. Error z value
> (Intercept) -6.8940 10.7802 -0.640
> log(dose) 0.9333 1.3441 0.694
> Pr(>|z|)
> (Intercept) 0.522
> log(dose) 0.487
>
> (Dispersion parameter for binomial family taken to be 1)
>
> Null deviance: 0.513878 on 3 degrees of freedom
> Residual deviance: 0.010356 on 2 degrees of freedom
> AIC: 6.5458
>
> Number of Fisher Scoring iterations: 5
>
> And the LD50 and CI transformed:
>
> print(EAUS.Aa)
> LD SE LCL UCL
> p = 0.5: 1614.444 3.207876 164.3822 15855.91
> p = 0.9: 6373.473 3.764879 474.1600 85669.72
>
> These are the values I get on SPSS (just replacing the values on R output) :
>
> Coefficients:
> Estimate Std. Error z value
> (Intercept) -6.8940 1.082 -6.373
> (dose) 2.149 0.311 6.918
>
> And the LD50 and CI transformed:
>
> LD LCL UCL
> p = 0.5: 1614.444 1198.932 1953.120
> p = 0.9: 6373.473 5145.767 9013.354
>
> So, please if somebody can help me with this, I'd be grateful. If working with those functions won't do it, I'll use another, the one you recommend.
>
> Thank you very much!
>
>
> Best wishes,
>
> Bianca
>
>
>
> PD. I've already googled it but there's no satisfactory answer.
>
>
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.