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