1
$\begingroup$

I want to report the results of a binomial logistic regression where I want to assess difference between the 3 levels of a factor (called System) on the dependent variable (called Response) taking two values, 0 and 1. My goal is to understand if the effect of the 3 systems (A,B,C) in System affect differently Response in a significant way. I am basing my analysis on this URL: https://stats.idre.ucla.edu/r/dae/logit-regression/

This is the result of my analysis:

> fit <- glm(Response ~ System, data = scrd, family = "binomial")
> summary(fit)

Call:
glm(formula = Response ~ System, family = "binomial", data = scrd)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8840   0.1775   0.2712   0.2712   0.5008  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)    3.2844     0.2825  11.626  < 2e-16 ***
SystemB  -1.2715     0.3379  -3.763 0.000168 ***
SystemC    0.8588     0.4990   1.721 0.085266 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.26  on 1023  degrees of freedom
Residual deviance: 376.76  on 1021  degrees of freedom
AIC: 382.76

Number of Fisher Scoring iterations: 6

Following this analysis I perform the wald test in order to understand whether there is an overall effect of System:

library(aod)

> wald.test(b = coef(fit), Sigma = vcov(fit), Terms = 1:3)
Wald test:
----------

Chi-squared test:
X2 = 354.6, df = 3, P(> X2) = 0.0

The chi-squared test statistic of 354.6, with 3 degrees of freedom is associated with a p-value < 0.001 indicating that the overall effect of System is statistically significant.

Now I check whether there are differences between the coefficients using again the wald test: 

# Here difference between system B and C:

> l <- cbind(0, 1, -1)
> wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 22.3, df = 1, P(> X2) = 2.3e-06



# Here difference between system A and C:

> l <- cbind(1, 0, -1)
> wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 12.0, df = 1, P(> X2) = 0.00052



# Here difference between system A and B:

> l <- cbind(1, -1, 0)
> wald.test(b = coef(fit), Sigma = vcov(fit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 58.7, df = 1, P(> X2) = 1.8e-14

My understanding is that from this analysis I can state that the three systems lead to a significantly different Response. Am I right? If so, how should I report the results of this analysis? What is the correct way?

|cite|improve this question
$\endgroup$
  • $\begingroup$ Interesting on your output: A vs C: glm part say p value = 0.085266, Wald's part say p = 0.00052. Something for A vs B. I am not expert on R, but this big difference is not explainable. $\endgroup$ – user158565 Nov 12 '18 at 1:28
  • $\begingroup$ I see more people do a likelihood ratio test instead of a wad test. Is there any reason you prefer the Wald test? $\endgroup$ – Demetri Pananos Nov 12 '18 at 1:28
  • $\begingroup$ Given the sample size > 1000, different tests should produce the very similar p-values. $\endgroup$ – user158565 Nov 12 '18 at 2:35
  • $\begingroup$ Thanks. The sample size is 1024 samples. I used the wald test as it is suggested here: stats.idre.ucla.edu/r/dae/logit-regression How can I perform the likelihood test? Also, is using contrasts a good method? If so what would it be the right R code to use for my case? $\endgroup$ – L_T Nov 12 '18 at 10:50
  • $\begingroup$ Could you put that 6 number in your question? # of positive and negative is each system. $\endgroup$ – user158565 Nov 12 '18 at 16:33
0
$\begingroup$

I think I managed to understand the procedure of using contrasts. Here you go my analysis:

> fit1 <- glm(Response ~ System, data = scrd, family = "binomial", contrasts = list(System=contr.treatment(levels(scrd$System), base=1)))
> summary(fit1)

Call:
glm(formula = Response ~ System, family = "binomial", data = scrd, 
    contrasts = list(System = contr.treatment(levels(scrd$System), 
        base = 1)))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8840   0.1775   0.2712   0.2712   0.5008  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)    3.2844     0.2825  11.626  < 2e-16 ***
SystemB  -1.2715     0.3379  -3.763 0.000168 ***
SystemC    0.8588     0.4990   1.721 0.085266 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.26  on 1023  degrees of freedom
Residual deviance: 376.76  on 1021  degrees of freedom
AIC: 382.76

Number of Fisher Scoring iterations: 6


> 
> 
>  
> fit2 <- glm(Response ~ System, data = scrd, family = "binomial", contrasts = list(System=contr.treatment(levels(scrd$System), base=2))) 
> summary(fit2)

Call:
glm(formula = Response ~ System, family = "binomial", data = scrd, 
    contrasts = list(System = contr.treatment(levels(scrd$System), 
        base = 2)))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8840   0.1775   0.2712   0.2712   0.5008  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)     2.0129     0.1853  10.860  < 2e-16 ***
SystemA   1.2715     0.3379   3.763 0.000168 ***
SystemC     2.1303     0.4512   4.722 2.34e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.26  on 1023  degrees of freedom
Residual deviance: 376.76  on 1021  degrees of freedom
AIC: 382.76

Number of Fisher Scoring iterations: 6

> 
> 
> 
> fit3 <- glm(Response ~ System, data = scrd, family = "binomial", contrasts = list(System=contr.treatment(levels(scrd$System), base=3))) 
> summary(fit3) 

Call:
glm(formula = Response ~ System, family = "binomial", data = scrd, 
    contrasts = list(System = contr.treatment(levels(scrd$System), 
        base = 3)))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8840   0.1775   0.2712   0.2712   0.5008  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)     4.1431     0.4113  10.072  < 2e-16 ***
SystemA  -0.8588     0.4990  -1.721   0.0853 .  
SystemB   -2.1303     0.4512  -4.722 2.34e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 411.26  on 1023  degrees of freedom
Residual deviance: 376.76  on 1021  degrees of freedom
AIC: 382.76

Number of Fisher Scoring iterations: 6

I would report the results as follows:

"Two binomial logistic regressions were performed using two different contrasts to ascertain the effects of the system type on the likelihood that participants report correct identifications. Results showed that participants’ accuracy was significantly lower for the system B compared to both system C (p < 0.001) and system A (p < 0.001), as well as that system C led to significantly higher identification accuracies than the system B (p < 0.001).”

Are these sentences correct? Also, is it ok to report only the p-values? Maybe the values for the coefficient beta?

Any comment?

|cite|improve this answer
$\endgroup$
  • 1
    $\begingroup$ It'd be better to post the questions as questions rather than as the final paragraph of an answer - please search the site first to see if they're already answered. $\endgroup$ – Scortchi - Reinstate Monica♦ Nov 10 '19 at 8:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.