1
$\begingroup$

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))

enter image description here

"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))

enter image description here

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.