2
$\begingroup$

I came across this quote from a paper by Pool and Doris (2021): "The 18–24 age group reported the highest levels of food insecurity, with 27.6% reporting being food insecure, as opposed to 8.2% of over 65-year-olds (odds ratio (OR) 2.45, 95% confidence interval (CI): 1.05; 5.75)." I'd like to understand how they calculated the probability of over-65s which was their reference group in the logistic regression they ran.

Example logistic regression output:

  Coefficients:
            Estimate Std.  Error z value Pr(>|z|)    
(Intercept) -2.457415 0.221487 -11.095  < 2e-16 ***
wave5       -0.510565 0.089884  -5.680 1.34e-08 ***
wave7       -0.360667 0.089068  -4.049 5.14e-05 ***
age_r2      -0.167394 0.163767  -1.022 0.306713    
age_r3      -0.664402 0.171821  -3.867 0.000110 ***
age_r4      -0.898270 0.161779  -5.552 2.82e-08 ***
age_r5      -1.704232 0.162821 -10.467  < 2e-16 ***
age_r6      -2.909736 0.183640 -15.845  < 2e-16 ***

I want to find out what the probability of the outcome variable is for the reference categories e.g. wave4 and age_r1. Apologies if this is straightforward but I've really struggled to find an answer. How can I do this in R?

Reference: Pool, U. and Dooris, M., 2021. Prevalence of food security in the UK measured by the Food Insecurity Experience Scale. Journal of Public Health,.

$\endgroup$
1
  • $\begingroup$ You could use the emmeans package to get the logit of each group. And then transform the logit to probability. $\endgroup$ May 13 at 1:06
1
$\begingroup$

From what you say, I assume that the combination (wave4, age_r1) means that all other $X$ variable values are 0. (Although it seems odd that there are only categories wave4, wave5 and wave7).

So, the estimated logit of the probability is $$-2.457415 +(0\times -0.510565) + (0 \times -0.360667) + (0 \times -0.167394) + (0 \times -0.664402) + (0 \times -0.898270) + (0 \times -1.704232) + (0 \times -2.909736)$$ $$= -2.457415 $$

Thus the estimated probability itself is $$\frac{\exp(-2.457415)}{1+\exp(-2.457415)}$$ $$= 0.07889799$$

That is the easiest one. You can use the same trick to get the probabilities for all other combinations too, although there will be some "$1\times$" multipliers.

You can see all of this using R the glm fitted object (call it "fit") by cbinding fit$fitted.values to your data set and scrolling through the result.

$\endgroup$

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.