2
$\begingroup$

I'm fairly new here so my apologies if I'm asking something obvious. My problem is the following:

I have a dataset in which I want to examine the interaction of a risk-factor with an intervention to predict a given outcome. The model looks as follows: 

summary(glm(Outcome~Risk.factor*Trial.arm,data=mydata,family="binomial"))

Which gives:

Estimate                        Estimate  Std. Error z value  Pr(>|z|)    

(Intercept)                      -1.5041     0.3191  -4.713 2.44e-06 ***
Trial.armTreated                  0.5773     0.4187   1.379  0.16796    
Risk.factorYes                    1.6582     0.6414   2.585  0.00973 ** 
Trial.armTreated:Risk.factorYes  -2.5232     1.0335  -2.441  0.01463 *  `

However, when I change my factor reference level (e.g. from using "Untreated" as reference, to using "Treated"), I get a sign reversal in the interaction and treatment term, as expected, but a completely different intercept and coefficient for the risk factor term. 

summary(glm(Outcome~Risk.factor*Trial.arm.rev,data=mydata,family="binomial"))

gives:

                                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)                          -0.9268     0.2710  -3.419 0.000628 ***
Trial.arm.revPlacebo                 -0.5773     0.4187  -1.379 0.167956    
Risk.factorYes                       -0.8650     0.8104  -1.067 0.285823    
Trial.arm.revPlacebo:Risk.factorYes   2.5232     1.0335   2.441 0.014631 *

I've tried tried the same with recoding all variables to binary, but the same happens. Does anyone know why this happens? The interpretation of the results seems quite different (it would appear to me).

Full Data here:

structure(list(Trial.arm = structure(c(1L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 
2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 
2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 
2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 
1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 
2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 
1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 
1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("Placebo", 
"Treated"), class = "factor"), Risk.factor = structure(c(2L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 
1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L), .Label = c("No", 
"Yes"), class = "factor"), Outcome = structure(c(1L, 1L, 2L, 
2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 
1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 
1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("No", 
"Yes"), class = "factor"), Risk.factor.binary = structure(c(1, 
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 
0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), .Label = c("No", "Yes")), 
    Outcome.binary = structure(c(0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 
    1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 
    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 
    0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
    0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 
    1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 
    1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0), .Label = c("No", 
    "Yes")), Trial.arm.binary = structure(c(0, 1, 1, 1, 1, 0, 
    0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 
    0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 
    1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 
    1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 
    1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 
    1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 
    1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 
    1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 
    1, 0), .Label = c("Placebo", "Treated")), Trial.arm.rev = structure(c(2L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
    1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 
    2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 
    1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 
    2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 
    2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 
    1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 
    1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 
    2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 
    2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
    1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("Treated", 
    "Placebo"), class = "factor"), Trial.arm.rev.bin = structure(c(1, 
    0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 
    1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 
    0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 
    1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 
    1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 
    0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 
    0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 
    1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 
    0, 0, 1, 0, 1, 0, 1), .Label = c("Treated", "Placebo"))), .Names = c("Trial.arm", 
"Risk.factor", "Outcome", "Risk.factor.binary", "Outcome.binary", 
"Trial.arm.binary", "Trial.arm.rev", "Trial.arm.rev.bin"), row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", 
"25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", 
"36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", 
"47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", 
"58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", 
"69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", 
"80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", 
"91", "92", "93", "94", "95", "96", "97", "98", "99", "100", 
"101", "102", "103", "104", "105", "106", "107", "108", "109", 
"110", "111", "112", "113", "114", "115", "116", "117", "118", 
"119", "120", "121", "122", "123", "124", "125", "126", "127", 
"128", "129", "130", "131", "132", "133", "134", "135", "136", 
"137", "138", "139", "140", "141", "142", "143", "144", "145", 
"146", "147", "148", "149", "150", "151", "152", "153", "154", 
"155", "156", "157", "158", "159", "160"), class = "data.frame") -> mydata
Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ The intercept is now for the Risk.FactorNo & Trial.armTreated combination. Before it was for the Risk.FactorYes & Trial.armPlacebo. So the intercept and the Risk.Factor terms do change to reflect this. The coefficient for Risk.FactorYes is now with reference to Risk.FactorNo & Trial.armTreated. $\endgroup$
    – Gavin Simpson
    Commented Jun 17, 2015 at 16:46
  • $\begingroup$ Hi Gavin, that of course makes sense and I see now why the intercept would be numerically different. I'm just really worried about the 'RiskFactorYes' increment (I only recoded trial.arm, not Risk.factor) relative magnitude. How come it's nearly twice as large and now significant? And by extension: how do you guard against that? For this case it's clear what the reference should be. But say it was male v female; setting the ref. category is arbitrary in that case. One could lead to a significant model parameter, the other not. Should you always check both? $\endgroup$
    – Hgremmels
    Commented Jun 17, 2015 at 17:29
  • $\begingroup$ Because you aren't comparing like with like. The intercept represents a certain combination of the two factors. The coef for RiskFactorYes is now the adjustment you make to the mean of the RiskFactorNo & Trial.armTreated group. $\endgroup$
    – Gavin Simpson
    Commented Jun 17, 2015 at 17:33
  • $\begingroup$ Thanks for the quick response :-). I still have difficulties wrapping my head around it though. Let me think out loud: The base odds for placebo|RiskfactorNo are exp(-1.5) = 0.22 and for treated|RiskFactorNo are exp(-0.92)=0.40. In both models the intercept for treatment is exp (0.57) which is the ratio of the two base odds in the models, that makes all perfect sense. Now in the placebo|RiskfactorNo reference model, the RiskFactorYes increases odds of outcome by exp(1.66)=5.26. In the treated|RiskfactorNo base model the effect of RiskFactorYes is exp(-0.87)=0.42, so an odds reduction. $\endgroup$
    – Hgremmels
    Commented Jun 17, 2015 at 18:08
  • $\begingroup$ But then the increment of Trial.armTreated:Risk.factorYes in the first model and of Trial.arm.revPlacebo:Risk.factorYes is the second is exactly the same, but in opposite sign; namely an odds ratio of exp(2.52)=12.4 on what seem completely different cumulative odds (up till then in model) $\endgroup$
    – Hgremmels
    Commented Jun 17, 2015 at 18:08

1 Answer 1

Reset to default
2
$\begingroup$

The predictions for each combination of Treated and Risk.factor are exactly the same under the two models. What I was trying to get at in the comments is that it is the parameterization of the model that is different; what is relative to what.

Assuming your first model is in m1 and your second model is in m2 and we have the following prediction data

pdat <- data.frame(Trial.arm = rep(c("Placebo", "Treated"), each = 2), 
                   Risk.factor = rep(c("No","Yes"), times = 2))
pdat <- transform(pdat, Trial.arm.rev = Trial.arm)

Then we get exactly the same odds for each combination of Trial.arm and Risk.factor with the two models:

> exp(predict(m1, pdat))
        1         2         3         4 
0.2222222 1.1666667 0.3958333 0.1666667 
> exp(predict(m2, pdat))
        1         2         3         4 
0.2222222 1.1666667 0.3958333 0.1666667

We can visualise the effects in each model using the effects package

library("effects")
plot(allEffects(m1))
plot(allEffects(m2))

Apart from the ordering, these are equivalent:

enter image description here

enter image description here

One reason for this is the dummy coding used by default in R. These "treatment" contrasts are not orthogonal but also they place the emphasis on a specific level or combination of levels of factors. If we switch to sum-to-zero contrasts, then the parameterization compares a level to the average of the other levels of that factor:

op <- options(contrasts = c("contr.sum", "contr.poly"))
m1 <- glm(Outcome ~ Risk.factor*Trial.arm,     data=mydata, family="binomial")
m2 <- glm(Outcome ~ Risk.factor*Trial.arm.rev, data=mydata, family="binomial")
rbind(M1 = unname(coef(m1)), M2 = unname(coef(m2)))

We note now that for both the original and the reversed model the coefficients have the same absolute value, the sign just changes.

> rbind(M1 = unname(coef(m1)), M2 = unname(coef(m2)))
        [,1]       [,2]       [,3]       [,4]
M1 -1.017112 -0.1983077  0.3421487 -0.6308064
M2 -1.017112 -0.1983077 -0.3421487  0.6308064

To some extent, your problem also comes from over-interpreting the main effects of the factors in the presence of an interaction term in which both participate. The "effect" of Risk.factor on the response depends entirely on what level of Treat.arm we are talking about. There is no single "effect" of Risk.factor == Yes in the model, but two effects:

  1. Risk.factor == Yes when Treat.arm == Placebo, and
  2. Risk.factor == Yes when Treat.arm == Treated

The principle of marginality says we should focus on the higher-order term, which makes it difficult (impossible?) to get one's head around what a main effect of a term that is marginal to a higher order interaction.

Improve this answer
$\endgroup$
0

Your Answer

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

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