thr3ads.net - R help - [R] Differences between SPSS and R on probit analysis [Feb 2017]

If this information is useful, please help other people find it:
Share via:

Biank M

2017-Feb-23 20:26 UTC

[R] Differences between SPSS and R on probit analysis

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]]

William Dunlap

2017-Feb-24 20:01 UTC

head link

[R] Differences between SPSS and R on probit analysis

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.

William Dunlap

2017-Feb-24 20:29 UTC

head link

[R] Differences between SPSS and R on probit analysis

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.

R help - Feb 2017 - Differences between SPSS and R on probit analysis

[R] Differences between SPSS and R on probit analysis

[R] Differences between SPSS and R on probit analysis

[R] Differences between SPSS and R on probit analysis