Here is some data to play with
set.seed(1e4)
color = factor(c(rep("blue",20), rep("red",20), rep("yellow",150)), levels=c("red","blue","yellow"))
x1 = rnorm(20+20+150,12,3)
x2 = c(rnorm(20,5), rnorm(20,5), rnorm(150,7))
x3 = c(rnorm(20,7), rnorm(20,5), rnorm(150,5))
x4 = c(rnorm(20,5), rnorm(20,7), rnorm(150,5))
RedOrNot = ifelse(color=="red","Red","NotRed")
df = data.frame(
color = color,
RedOrNot = RedOrNot,
x1 = x1,
x2 = x2,
x3 = x3,
x4 = x4
)
We are interested in the color red in particular. We would like to understand which explanatory variables (which among x1,x2,x3 and x4) affects the probability red "differentially that it affects the probability of outcome blue or yellow". The question might be a little unclear but I hope the following will make my point very obvious.
Below graphs show 95% bootstrap confidence intervals.
p1 = ggplot(df,aes(x=color,y=x1)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
p2 = ggplot(df,aes(x=color,y=x2)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
p3 = ggplot(df,aes(x=color,y=x3)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
p4 = ggplot(df,aes(x=color,y=x4)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
multiplot(p1,p2,p3,p4,layout=matrix(1:4,ncol=2))
"red" seem to differ from the two other groups only in respect to x4. I would be tempted to group yellow and blue into a single group called NotRed as shown below.
p1 = ggplot(df,aes(x=RedOrNot,y=x1)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
p2 = ggplot(df,aes(x=RedOrNot,y=x2)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
p3 = ggplot(df,aes(x=RedOrNot,y=x3)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
p4 = ggplot(df,aes(x=RedOrNot,y=x4)) + stat_summary(fun.data=mean_cl_boot) + coord_flip()
multiplot(p1,p2,p3,p4,layout=matrix(1:4,ncol=2))
Because we sampled way more yellow than blue, it seems that Red differ from the other colors for x2 as well as x4 but this is misleading. I would rather not consider x2 as having an effect of interest when I see the first series of plots as I am looking for variables that make red different than both other groups.
Let's make a test (a categorical MCMCglmm) with color as response variable.
k = length(levels(df$color))
I = diag(k-1)
J = matrix(rep(1, (k-1)^2), c(k-1, k-1))
priors = list(R = list(fix=1, V=(1/k) * (I + J), n = k - 1))
m = MCMCglmm(color ~ -1 + trait:(x1) + trait:(x2) + trait:(x3) + trait:(x4) , rcov = ~ us(trait):units, data = df, family = "categorical", prior = priors, verbose = FALSE, burnin = 10000, nitt = 60000, thin = 50)
summary(m)
post.mean l-95% CI u-95% CI eff.samp pMCMC
traitcolor.blue:x1 -0.30261 -0.60335 0.02378 72.04 0.048 *
traitcolor.yellow:x1 -0.02728 -0.24525 0.19561 129.73 0.832
traitcolor.blue:x2 -0.60401 -1.20394 0.04258 32.72 0.118
traitcolor.yellow:x2 2.05721 1.35060 2.94990 26.03 <0.001 ***
traitcolor.blue:x3 2.27062 1.36630 3.27224 30.57 <0.001 ***
traitcolor.yellow:x3 -0.30359 -0.97326 0.27461 86.95 0.306
traitcolor.blue:x4 -1.28126 -2.11725 -0.55014 41.71 <0.001 ***
traitcolor.yellow:x4 -1.41233 -2.05106 -0.80523 68.91 <0.001 ***
Would it be wise to consider only the highest p.value for each explanatory variable to decide (relative to a threshold) whether a variable affects the probability of the outcome red differentially than both blue and yellow?
How can I go about testing what variables affect the probability of the outcome red "differentially than both blue and yellow"?
Similarly, when performing an LDA, I would have wished x4 to pop up as important driver of the first axis but this is not made obvious due to x2 again.
lda(data=df,RedOrNot~x1+x2+x3+x4)$scaling
LD1
x1 0.004390466
x2 -0.571945258
x3 -0.241371706
x4 0.784905799
lda(data=df,color~x1+x2+x3+x4)$scaling
LD1 LD2
x1 -0.0361047 0.03607743
x2 -0.8110918 -0.09303029
x3 0.4472190 -0.69882190
x4 0.3820449 0.73839111
