How to determine if two categorical variables are dependent while controlling for a 3rd categorical? I have 3 categorical variables: country, gender, and liked (whether the user liked the content or not). Using Chi-squared I see that 'liked' is dependent on country, that 'liked' is dependent on gender, and gender is dependent on country (ran chi-square separately on each pair of these variables and p was < 0.05). 
I want to check whether gender is a proxy to country - so essentially whether the graph is country->gender->liked or whether country impacts both the gender and has a direct effect on 'liked'. 
What test can I run to see if liked is independent of country when controlling for gender? (again the 3 variables are categorical) 
 A: This sounds like a logistic regression case to me with "liked" as outcome variable.
But you might also have a look at the cotaplot:
library(vcd)
head(dat <- as.data.frame(Titanic))
# 
# I. structured contingency table         # II. mosaic plot
# (tab <- xtabs(Freq ~ ., data=dat));       mosaic(tab, shade=T)                            # significant connection between any?
structable(~ Class + Survived, data=dat); mosaic(~ Class + Survived, data=dat, shade=T)   # significant connection between x & y?
structable(~ Sex   + Survived, data=dat); mosaic(~ Sex   + Survived, data=dat, shade=T)   # significant connection between x & y?
structable(~ Age   + Survived, data=dat); mosaic(~ Age   + Survived, data=dat, shade=T)   # significant connection between x & y?
# 
cotabplot(~ Survived + Sex | Class, data=dat, shade=T)  # significant connection between x & y in all z?
# 
# a) -> Independence can be clearly rejected (p < 2.22*10^-16). 
#       In 1st class surprisingly many people survived, 
#       while in 3rd class and among the crew surprisingly many people died.
# b) -> Independence can be clearly rejected (p < 2.22*10^-16). 
#       Women survived much more, so they were probably saved first.
# c) -> Independence can be clearly rejected (p < 4.7*10^-6)
#       Children survived more, so they were probably saved first.
# d) -> Independence can be clearly rejected in each class. Women were always favoured; 
#       it seems, that they were especially favoured in class 1, less so in class 3 or crew.

A: It's been a while but I try it anyway ...
My notes say concerning 2 categorical variables:
# Independence testing & comparison of samples
(tab.dat <- table(dat$A, dat$B))
assocplot(tab.dat)
(prop.dat <- prop.table(tab.dat,1) # B geg. A
(prop.dat <- prop.table(tab.dat,2) # A geg. B

# X^2 has more power in certain cases:
chisq.test(correct=F) # = 8.33
chisq.test(correct=T) # = 7.35, Yates correction at n<60

# Fisher only covers discrete significance levels, 
# but is more accurate:
fisher.test()  
                  
# For connected samples:
mcnemar.test()
binom.test() # more exact at ~ Bin(n12+n21,0.5)

And for logistic regression I found this example with residual analysis
(http://stat.ethz.ch/Teaching/Datasets/WBL/kevlar49.dat):
head(dat <- read.table("kevlar49.dat", header=T))
str(dat)
dat$Spool <- factor(dat$Spool)
dd <- dat[dat$Stress >= 24, ]
# 
a) Adapt a linear model and do a residual analysis:
# log(Failure) = ßo + ß_1Stress + ß_2Spool + E, E~N(0, sigma^2)
# In the following subtasks we want to do the Weibull regression.
require(regr0)
(r.dat <- regr(log(Failure) ~ Stress + Spool, data=dd))
plot(r.dat, sequence=T)
# All input variables are significant.
# residue analysis:
# . TA: the smoothing curve is slightly curved
# . Scatter diagram: the smoothing curve is falling
# . QQ-Plot: the QQ-Plot does not fit optimally, 
    # is rather short-tailed
# . Leverage plot: few data have a larger lever arm.
# . Input variable residuals (not shown in the sample solution): 
    # acceptable
# . Independence: is violated
# The assumptions are not fulfilled overall.

I don't remember exactly, how I got to the conclusion of violated independence by looking at the residual analysis plots ... maybe someone else can help out?
A: In more abstract, model independent terms, you're trying to show that 
$P(L|C,G)=P(L|G)$
(if you know the gender, then knowing the country doesn't give you any extra information)
You can train a classifier (could be logistic regression, but could also try random forest, xgboost, neural network, etc) which knows only about Gender, and another one which knows about Country and Gender. You can then, under cross-validation see whether the latter model outperforms the former. 
If your hunch is correct, you should find that they outperform each other about equally frequently. 
The main downfall of this method, is if $P(L|C,G)\simeq P(L|G)$, i.e. knowing country gives you a small amount of extra information over knowing only gender.
Then you might find, that the model which includes both outperforms the model which includes gender only, 80% of the time. Or it might outperform it every time, but with tiny differences in cross-entropy. Then it can be quite difficult to determine whether it's "worth" including country or whether its contribution is so marginal that it's not worth it. 
