9
$\begingroup$

I have following table of values: 

25  75
38  162

The odds ratio is 0.7037 and log(OR) is -0.3514. For a contingency table with values a, b, c, and d, the variance of log(OR) is given by

(1/a + 1/b + 1/c + 1/d)

How can I calculate the p.value of log(OR) from this data in R (whether it is significantly different from 0)?

Improve this question
$\endgroup$

4 Answers 4

Reset to default
11
$\begingroup$

You can use Fisher's exact test, which inputs a contingency table and outputs a p-value, with a null hypothesis that the odds ratio is 1 and an alternative hypothesis that the odds ratio is not equal to 1.

(tab <- matrix(c(38, 25, 162, 75), nrow=2))
#      [,1] [,2]
# [1,]   38  162
# [2,]   25   75
fisher.test(tab)
# 
#   Fisher's Exact Test for Count Data
# 
# data:  tab
# p-value = 0.2329
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  0.3827433 1.3116294
# sample estimates:
# odds ratio 
#  0.7045301

In this case the p value is 0.23.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for a clever way to determine the p value. Chi-square test can also be used in a similar manner. $\endgroup$
    – rnso
    Jun 13, 2015 at 17:09
  • $\begingroup$ @rnso for sure, though Fisher's exact test is preferred over Chi-square when you have small cell sizes in your contingency table. $\endgroup$
    – josliber
    Jun 13, 2015 at 17:16
  • 4
    $\begingroup$ This is a longstanding myth but is unfortunately not true. The ordinary Pearson $\chi^2$ provides more accurate $P$-values than Fisher's so-called "exact" test even when expected frequencies are as low as 1.0. $\endgroup$
    – Frank Harrell
    Jun 14, 2015 at 12:14
  • $\begingroup$ could you say a bit more about this @FrankHarrell? I know the $\chi^2$ would be an asymptotic result, whereas the Fisher's exact test relies on the exact distribution, how is the $p$-value more "accurate" using the asymptotic method? $\endgroup$
    – bdeonovic
    Jun 15, 2015 at 11:58
  • 1
    $\begingroup$ See extensive comments about this on the site. Briefly, the P-values from Fisher's test are too large. The mean absolute error in P-values from Pearson's test is smaller. Fisher's is only "exact" in the sense that P-values are "guaranteed" not to be too small. $\endgroup$
    – Frank Harrell
    Jun 16, 2015 at 12:23
10
$\begingroup$

Another way to do it (other than Fisher's exact test) is to put the values into a binomial GLM:

d <- data.frame(g=factor(1:2),
                s=c(25,75),
                f=c(38,162))
g <- glm(s/(s+f)~g,weights=s+f,data=d,
    family="binomial")
coef(summary(g))["g2",c("Estimate","Pr(>|z|)")]
##   Estimate   Pr(>|z|) 
## -0.3513979  0.2303337

To get the likelihood ratio test (slightly more accurate than the Wald $p$-value shown above), do

anova(g,test="Chisq")

which gives

##      Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL                     1     1.4178         
## g     1   1.4178         0     0.0000   0.2338

(LRT $p=0.2338 \approx$ Wald $p=0.2303337 \approx$ Fisher $p=0.2329$ in this case because the sample is fairly large)

Improve this answer
$\endgroup$
5
$\begingroup$

It's better to generalize the solution and use the likelihood ratio $\chi^2$ test from a statistical model such as the logistic model. The LR test provides fairly accurate $P$-values. This also handles cases where you need to test more than one parameter, e.g., 3-group problems, continuous effects that are nonlinear, etc. The LR test for the overall model (which is all that's needed in this example since there are no adjustment variables) may be easily obtained in base R or using the rms package, e.g.

f <- lrm(y ~ groups, weights=freqs)
f  # prints LR chi-sq, d.f., P, many other quantities

Here the nested models are this model and an intercept-only model.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I could find that LR test (lrtest) is used for comparing nested models. How can we use it here? Could you write a line of R code for it? $\endgroup$
    – rnso
    Jun 15, 2015 at 2:36
  • $\begingroup$ for what it's worth this is more or less the same statistical approach (although with a better explanation) as in my answer above. lrm() has different defaults, output formats, etc., but the statistical model (IIUC) is the same as glm(...,family="binomial") $\endgroup$
    – Ben Bolker
    Jun 24, 2015 at 21:51
2
$\begingroup$

An even simpler code is to use directly a poisson regression with

library(tidyverse)
library(broom)

d <- data.frame(g=factor(1:2),
                s=c(25,75),
                f=c(38,162))

d_long = pivot_longer(d, cols = 2:3)

tidy(glm(value ~ name*g, data=d_long, family = 'poisson'))

Some explanations

If you look at the poisson model, you see that it models the probability $\pi$ as

\begin{equation} \ln(\pi_{i}) = \mathbf{X_{i}} \beta \end{equation}

the log of probabilities is a simple linear combinations of predictors. This is the same equation as the logistic equation which is called "logit"

\begin{equation} \text{logit}(\pi_{i}) = \ln \left( \frac{\pi_{i}}{1 - \pi_{i}} \right) = \mathbf{X_{i}} \beta \end{equation}

You can test this equivalence with a simple model. You start with this linear model

\begin{equation} y_{i} = 2u_{i} + \epsilon_{i} \end{equation}

set.seed(123)
n = 100
u = rbinom(n, 1, 0.2)
y = u*2 + rnorm(n, 0, 0.5)

y_prob = plogis(y) # transform log(odds) into a proba
y_bin = rbinom(n, 1, prob = y_prob)

You can directly model the data with a logistic regression

glm(y_bin ~ u, family = "binomial") 

Coefficients:
(Intercept)            u  
   -0.09764      2.93085

This is the log of odds, which we can simply retrieve with a table

 data.frame(table(u, y_bin)) %>% mutate(p = Freq/sum(Freq))

$$ \begin{aligned}[ht] \begin{array}{rllrr} \hline & u & y\_bin & Freq & p \\ \hline 1 & 0 & 0 & 43 & 0.43 \\ 2 & 1 & 0 & 1 & 0.01 \\ 3 & 0 & 1 & 39 & 0.39 \\ 4 & 1 & 1 & 17 & 0.17 \\ \hline \end{array} \end{aligned} $$

The coefficient of your logistic regression is the log odds-ratio of this table

\begin{equation} \ln \left( \frac{\frac{0.17}{0.39}}{\frac{0.01}{0.43}} \right) = 2.930852 \end{equation}

log((0.17/0.39) / (0.01/0.43))
2.930852

This is exactly the same as modelling the table directly with a poisson regression

Just prepare the table

df_poiss = data.frame(table(u, y_bin))

Run the poisson, and the interaction will give you the log odds-ratio between the outcome and the predictor u

glm(Freq ~ y_bin*u, data = df_poiss, family = "poisson")

Coefficients:
(Intercept)       y_bin1           u1    y_bin1:u1  
    3.76120     -0.09764     -3.76120      2.93085

The problem is that you can always transform dummy variables into counts for the purpose of using a poisson model. But you often cannot use aggregated counts to revert back to dummies for using logistic regression.

A final note if you look closely at the two equations, the logistic model (the logit transformation) once exponentiate will give you odds, while when exponentiating the poisson, will give you a probability (or risk). However, when you take an interaction with the poisson, then you have a ratio of odds, which is the same as the logistic model.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Some explanation on the use of poisson vs binomial will be helpful. $\endgroup$
    – rnso
    Feb 27 at 18:07
  • $\begingroup$ Very well explained. $\endgroup$
    – rnso
    Feb 28 at 12:48
  • $\begingroup$ thanks no worries. Hope it helps $\endgroup$
    – giac
    Feb 28 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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