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_d...@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@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. ______________________________________________ R-help@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.