[R] Differences between SPSS and R on probit analysis

[William Dunlap]

Another model specification equivalent to
    cbind(afflicted, total-afflicted) ~ ...
is the ratio you had accompanied by the total as the 'weights' argument
    afflicted/total ~ ..., weights=total
Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Fri, Feb 24, 2017 at 12:01 PM, William Dunlap <wdunlap at tibco.com> wrote:
> 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.