Logistic regression model comparison

I want to use logistic regression to look at how a continuous variable that I have measured for a sample of participants (how long they looked at a product) is impacting a dichotomous variable (whether they bought it or not, Y or N). There are several products I am interested looking at, which means I will have to run separate LR's for each product, I believe.

I am mainly interested in comparing the relative impact that the IV has on purchases of each product. I thought that the most intuitive way to do this would be using the Odds Ratios. I have two questions about this.

1) Is it appropriate to compare the odds ratios between the models?

2) Will the model's predictive ability/performance be important to interpret if I am only interested in the odds ratios? What does a poorly performing model mean for the interpretation of the odds ratio?

Thanks for any help, it is very much appreciated.

• You have multiple outcomes for each person, so I would use a multilevel model to do this. You then test an interaction effect across products. Note that interactions with logistic models are a bit tricky though. Mar 26 '13 at 17:43
• Thanks for the response. Are you referring to conditional or multinomial LR? Participants were not constrained to choose only one product, so I don't think this would work. Mar 26 '13 at 18:52
• Ah, perhaps I misunderstood. Were all products presented simultaneously? Could they choose more than one, or just one? Mar 26 '13 at 22:13
• Yes, products were presented simultaneously and more than one could be chosen. In fact, multiple products were chosen by a majority of the participants. So I think different logistic models for each product is the only option. I am just unsure if comparing these models' odds ratios is appropriate, and whether model performance is important if I am only interested in said odds ratios. Mar 27 '13 at 13:53
• I think multilevel binary logistic regression is the way to go. This will estimate an OR for each product associated with gaze length, and let you test the difference between the ORs. I'll respond at more length later, if someone else doesn't. Mar 27 '13 at 14:44

Here's how I'd do it.

First, I'll make up some data

require(lme4)
set.seed(1234)
time = rnorm(500) + 10
id <- rep(c(1:100), 5)

product <- rep(1:5, 100)
product <- product[order(product)]
selected <- rbinom(prob=0.3, size =1, n=500)

d <- as.data.frame(cbind(time, id, product, selected))
d$product <- as.factor(d$product)
d$id <- as.factor(d$id)


A couple of logistic regressions, to check things out, for individual products.

logistic1 <- glm(selected ~ time, data=subset(d, product==1), family=binomial)
summary(logistic1)

logistic2 <- glm(selected ~ time, data=subset(d, product==2), family=binomial)
summary(logistic2)

> Call: glm(formula = selected ~ time, family = binomial, data =
> subset(d,
>     product == 1))
>
> Deviance Residuals:     Min      1Q  Median      3Q     Max
> -1.173  -1.092  -1.054   1.266   1.338
>
> Coefficients:
>             Estimate Std. Error z value Pr(>|z|) (Intercept) -1.02762    1.99176  -0.516    0.606 time         0.08398    0.20116   0.417    0.676
>
> (Dispersion parameter for binomial family taken to be 1)
>
>     Null deviance: 137.63  on 99  degrees of freedom Residual deviance: 137.45  on 98  degrees of freedom AIC: 141.45
>
> Number of Fisher Scoring iterations: 3
>
> >
> > logistic2 <- glm(selected ~ time, data=subset(d, product==2), family=binomial)
> > summary(logistic2)
>
> Call: glm(formula = selected ~ time, family = binomial, data =
> subset(d,
>     product == 2))
>
> Deviance Residuals:
>     Min       1Q   Median       3Q      Max
> -1.0646  -1.0114  -0.9946   1.3446   1.4115
>
> Coefficients:
>             Estimate Std. Error z value Pr(>|z|) (Intercept) -0.85492    2.00998  -0.425    0.671 time         0.04474    0.19895   0.225    0.822
>
> (Dispersion parameter for binomial family taken to be 1)
>
>     Null deviance: 134.60  on 99  degrees of freedom Residual deviance: 134.55  on 98  degrees of freedom AIC: 138.55
>
> Number of Fisher Scoring iterations: 4


Then two lmer models, one with and one without the product * time interaction term.

lmerNoInt <- lmer(selected ~ time + product  + (1 | id)  , data=d, family=binomial)
summary(lmerNoInt)

lmerWithInt <- lmer(selected ~ time + product  + time * product + (1 | id)  , data=d,     family=binomial)
summary(lmerwithInt)


Here's the output:

> lmerNoInt <- lmer(selected ~ time + product  + (1 | id)  , data=d, family=binomial)
> summary(lmerNoInt)
Generalized linear mixed model fit by the Laplace approximation
Formula: selected ~ time + product + (1 | id)
Data: d
AIC   BIC logLik deviance
698.1 727.6 -342.1    684.1
Random effects:
Groups Name        Variance Std.Dev.
id     (Intercept)  0        0
Number of obs: 500, groups: id, 100

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.4875     0.8911  -1.669   0.0951 .
time          0.1307     0.0881   1.483   0.1380
product2     -0.2317     0.2877  -0.805   0.4206
product3      0.4028     0.2863   1.407   0.1594
product4      0.1015     0.2846   0.357   0.7213
product5      0.2642     0.2845   0.929   0.3531
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) time   prdct2 prdct3 prdct4
time     -0.974
product2 -0.094 -0.066
product3 -0.073 -0.088  0.498
product4 -0.116 -0.045  0.499  0.502
product5 -0.125 -0.036  0.498  0.501  0.503
>
>
> lmerWithInt <- lmer(selected ~ time + product  + time * product + (1 | id)  , data=d, family=binomial)
> summary(lmerwithInt)
Generalized linear mixed model fit by the Laplace approximation
Formula: selected ~ time + product + time * product + (1 | id)
Data: d
AIC BIC logLik deviance
701.6 748 -339.8    679.6
Random effects:
Groups Name        Variance   Std.Dev.
id     (Intercept) 1.6416e-10 1.2812e-05
Number of obs: 500, groups: id, 100

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.02771    1.99178  -0.516    0.606
time           0.08399    0.20117   0.418    0.676
product2       0.17265    2.82970   0.061    0.951
product3      -3.96985    3.09480  -1.283    0.200
product4       1.79796    2.77276   0.648    0.517
product5      -0.59914    2.70511  -0.222    0.825
time:product2 -0.03923    0.28293  -0.139    0.890
time:product3  0.43306    0.30804   1.406    0.160
time:product4 -0.16911    0.27813  -0.608    0.543
time:product5  0.08716    0.27163   0.321    0.748

Correlation of Fixed Effects:
(Intr) time   prdct2 prdct3 prdct4 prdct5 tm:pr2
time        -0.995
product2    -0.704  0.700
product3    -0.644  0.640  0.453
product4    -0.718  0.715  0.506  0.462
product5    -0.736  0.733  0.518  0.474  0.529
time:prdct2  0.707 -0.711 -0.995 -0.455 -0.508 -0.521
time:prdct3  0.650 -0.653 -0.457 -0.996 -0.467 -0.478  0.464
time:prdct4  0.720 -0.723 -0.506 -0.463 -0.995 -0.530  0.514
time:prdct5  0.737 -0.741 -0.519 -0.474 -0.529 -0.994  0.527
tm:pr3 tm:pr4
time
product2
product3
product4
product5
time:prdct2
time:prdct3
time:prdct4  0.472
time:prdct5  0.484  0.536
>


The logistic regression told us that the coefficients for product 1 and product 2 were 0.08398 and 0.04474. In the lmer model with interaction, the effects were 0.08399 for time (which is for product 1) and for time * product 2, -0.03923. This is the difference between products 1 and 2. 0.08398 - 0.03923 = 0.04475, which is the estimate of the effect of time for product 2. The so the standard error on each interaction term is the difference between the effect of time for products 1 and 2. Same for all the others.

Then use the anova() function to get an omnibus test of the difference between the models. That tells you if the betas differ across products (and the betas are the log(OR)s.)

anova(lmerNoInt, lmerWithInt)


That's the test of whether the odds ratios of all 5 products differ.

> anova(lmerNoInt, lmerWithInt)
Data: d
Models:
lmerNoInt: selected ~ time + product + (1 | id)
lmerWithInt: selected ~ time + product + time * product + (1 | id)
Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
lmerNoInt    7 698.14 727.64 -342.07
lmerWithInt 11 701.61 747.97 -339.81 4.5268      4     0.3394


You could also do this with Lavaan, and get the same result, by setting up a multivariate logistic regression and adding a constraint. Edit: No you can't. You can do it with the (non-free) Mplus and you can do it with openMx, if you're cleverer than me.