# Correlation between two binary variables within one categorical variable

## The Problem:

I have measured two binary variables within 1 categorical variable with 5 levels.

Initially, I thought I'd be able to use Fisher's Exact test or some $$N \times M \times K$$ version of it. However I have only found references where the dimensions of the table are specifically $$2 \times 2 \times 2$$ (see the "hypergea" package in R as an example). In my case it is a $$5 \times 2 \times 2$$.

I was wondering if there are other approaches that would be appropriate.

Here's an example of what I'm working with in R.

set.seed(66)
group1<-as.data.frame(rbind(cbind(rep("GROUP_1",times=50),sample(x=c(1:50),size=50,replace=TRUE),rbinom(50,1,0.25),rep(NA,times=50)),cbind(rep("GROUP_1",times=50),sample(x=c(51:100),size=50,replace=TRUE),rep(NA,times=50),rbinom(50,1,0.75))))
group2<-as.data.frame(rbind(cbind(rep("GROUP_2",times=50),sample(x=c(101:150),size=50,replace=TRUE),rbinom(50,1,0.05),rep(NA,times=50)),cbind(rep("GROUP_2",times=50),sample(x=c(151:200),size=50,replace=TRUE),rep(NA,times=50),rbinom(50,1,0.95))))
group3<-as.data.frame(rbind(cbind(rep("GROUP_3",times=50),sample(x=c(201:250),size=50,replace=TRUE),rbinom(50,1,0.50),rep(NA,times=50)),cbind(rep("GROUP_3",times=50),sample(x=c(251:300),size=50,replace=TRUE),rep(NA,times=50),rbinom(50,1,0.50))))
group4<-as.data.frame(rbind(cbind(rep("GROUP_4",times=50),sample(x=c(301:350),size=50,replace=TRUE),rbinom(50,1,0.67),rep(NA,times=50)),cbind(rep("GROUP_4",times=50),sample(x=c(351:400),size=50,replace=TRUE),rep(NA,times=50),rbinom(50,1,0.33))))
group5<-as.data.frame(rbind(cbind(rep("GROUP_5",times=50),sample(x=c(301:350),size=50,replace=TRUE),rbinom(50,1,0.20),rep(NA,times=50)),cbind(rep("GROUP_5",times=50),sample(x=c(351:400),size=50,replace=TRUE),rep(NA,times=50),rbinom(50,1,0.20))))

testdata<-rbind(group1,group2,group3,group4)
names(testdata)<-c("GROUP","ID","Var1_","Var2_")

I'll note here that Var1_ and Var2_ were measured from different individuals, though the individuals were sampled from the same population. That's why I placed the NA's in the data frame. Also, note that some of the data may be pseudoreplicated by ID - I'll come back to that later.

I'm looking for a way of analyzing this problem.

Specifically I want to know:

1. Does the state of Var1_ significantly predict the state of Var2_? What is the relationship and how strong is it?

2. Is this relationship significantly different between GROUP levels? How different is it (effect size)?

## My approach so far:

Some searching indicated that a GLM might be appropriate. However, as far as I can tell, the GLM function requires the V1 and V2 variables to be observations from the same sample. So, just to see how it looks, I'll take the NA's out and run a GLM.

group1<-as.data.frame(cbind(rep("GROUP_1",times=50),rbinom(50,1,0.25),rbinom(50,1,0.75)))
group2<-as.data.frame(cbind(rep("GROUP_2",times=50),rbinom(50,1,0.05),rbinom(50,1,0.95)))
group3<-as.data.frame(cbind(rep("GROUP_3",times=50),rbinom(50,1,0.50),rbinom(50,1,0.50)))
group4<-as.data.frame(cbind(rep("GROUP_4",times=50),rbinom(50,1,0.67),rbinom(50,1,0.33)))
group5<-as.data.frame(cbind(rep("GROUP_5",times=50),rbinom(50,1,0.20),rbinom(50,1,0.20)))

testdata<-rbind(group1,group2,group3,group4)
names(testdata)<-c("GROUP","Var1_","Var2_")

glm_obj_1 <- glm(Var1_ ~ Var2_, data=testdata, family=binomial(link="logit"))

summary(glm_obj)

Ok... I think that might be what I'm after, but I'd like to know if the relationship is different between GROUPs.

glm_obj_2 <- glm(Var1_ ~ Var2_* GROUP, data=testdata, family=binomial(link="logit"))
summary(glm_obj_2)

Ok! Is this a reasonable approach? Removing the NA's seems questionable to me, but I'm not sure how else to deal with that in the GLM specification.

## Extra Credit:

I mentioned the pseudoreplication before. I assume this could be controlled for with a random effect? I'll try a GLMM implemented in LME4... This time I'll add fewer ID levels to increase the odds of pseudoreplication.

library(lme4)
group1<-as.data.frame(cbind(rep("GROUP_1",times=50),sample(x=c(1:5),size=50,replace=TRUE),rbinom(50,1,0.25),rbinom(50,1,0.75)))
group2<-as.data.frame(cbind(rep("GROUP_2",times=50),sample(x=c(6:10),size=50,replace=TRUE),rbinom(50,1,0.05),rbinom(50,1,0.95)))
group3<-as.data.frame(cbind(rep("GROUP_3",times=50),sample(x=c(11:15),size=50,replace=TRUE),rbinom(50,1,0.50),rbinom(50,1,0.50)))
group4<-as.data.frame(cbind(rep("GROUP_4",times=50),sample(x=c(16:20),size=50,replace=TRUE),rbinom(50,1,0.67),rbinom(50,1,0.33)))
group5<-as.data.frame(cbind(rep("GROUP_5",times=50),sample(x=c(21:25),size=50,replace=TRUE),rbinom(50,1,0.20),rbinom(50,1,0.20)))

testdata<-rbind(group1,group2,group3,group4)
names(testdata)<-c("GROUP","ID","Var1_","Var2_")

glmer_obj<-glmer(Var1_ ~ Var2_ * GROUP + (1|ID), data=testdata, family=binomial(link="logit"))

summary(glm_obj)

Well, it fails to converge, but for illustrative purposes, I hope it suffices. Would this be an appropriate way to tackle this type of data?

I'm eager to hear your responses, suggestions, criticisms, etc.

• How about Cochran–Mantel–Haenszel statistics ?
– MJW
Apr 19 '17 at 10:22
• Every ID is one unit (e..g, person), right? It seems like you don't have any units for which both Var1_ and Var2_ were observed: Then I don't think that you can say much about their relationship... Apr 19 '17 at 12:39
• @hplieninger That's what I was afraid of. Though, I don't understand why. A concrete example: You measure the prevalence of smoking (1= smoke, 0 = don't smoke) among 5 age groups. Each group has a different probability of smoking. Then you measure the prevalence of lung cancer (again, assume a binary response) in the same age groups, though in different individuals. Aren't they both estimating the value of a trait of each of those populations? Why must the data come from the same people if we know the people were sampled from the same population? Apr 23 '17 at 10:52
• Because the relationship at the group level may or may not be the same as the relationship at the individual level (which you are probably interested in). See ecological fallacy. Apr 24 '17 at 10:41
• I think it would help if your question was more concrete. What are these groups? experimental groups? natural groups? Why is either Var1_ OR Var2_ observed? Why are some IDs measured more than once and others not. Apr 24 '17 at 10:46