7
$\begingroup$

I'm trying to figure out what is the most appropriate statistic test for use with my data and was hoping for some advice.

The primary data consists of a binary independent variable (patient test - positive or negative; there are some missing values which are being excluded for now) for which we want to determine if it is predictive for a given clinical outcome.

The outcome is ordinal (negative or grade 1 - 4), however we are interested in any positive outcome so this could be simplified to a binary outcome (negative or positive), although if the test is only predictive for some grades this would be good to know.

In the simplified case, I could use a chi-square test. I'm less certain how to deal with the ordinal outcome.

Additionally, each patient has multiple tests and outcomes measured at different times, but the number of repeats and the times are not the same for all patients.

There are also a number of other data pertaining to each patient, some of which are binary (eg. sex, one of 2 research centers, certain other test results), ordinal (different test results) and continuous (age, yet some more test results which - in some cases complicated by a continuous result or N/A).

If the outcome were continuous I think I'd be looking at using a mixed effects linear model so I could determine if there were any effects or interactions with the various secondary independent variables. Is there something equivalent I can do for ordinal and/or binary outcome?

Thanks for your advice! Ric

Edit1 (in response to answer below by Macro): I did read your linked answer for more discussion re: random effects vs GEE. However, I'm honestly still unclear on when you would want an individual vs a population based coefficient.

To address my data specifically, we want to determine how well, if at all, the hypothesized predictor can predict the outcome. If we find it is a good predictor then we would want to be able to say when testing new patients in the future that given they are positive for test X, they are Z times more likely to become positive for outcome Y -- this will affect the course of treatment.

So would I want a random effects or a GEE model to answer this question? What would be a similar scenario where I would want the other model?

As a preliminary experiment I did a Fisher's Exact test on X vs Y (ignoring the fact that many observations are repeated measures) and found evidence for a possible trend but not a statistically significant effect (p ~ 0.1).

I repeated this after splitting the data into 2 sets based on research centre. Centre 1: p=1, odds ratio ~ 1. Centre 2: p=0.02, odds ratio ~ 3. Other independent factors may influence whether X predicts Y or may correlate with Centre.

Is logistic regression still an appropriate method to explorer these effects on the predictive effect of X? If it just tells me that X is not predictive but some other factor is, while interesting, it would not answer our primary question. If it tells us that X is predictive, but only for Centre 2, and especially for older patients at early times since treatment, that is the sort of thing we are looking for.

Thanks!

Edit2 (interpreting results?): So I've been playing with this the last few days but still having trouble making sense of it. I've found some helpful texts on logistic regression using glm() but nothing on lme4 that I was able to fully grok. Perhaps if you know of a good tutorial on the subject...

In any case, to fill in some details on the data set here's a table summarizing the variable structures:

variable  type      levels   missing values
id        factor      76           0
inst      factor       2           0
d1        Date        NA           0
d2        integer     NA           0
d3        Date        NA           0
d4        integer     NA           0
a         integer     NA           2
b         integer     NA           6
c         factor       2           2
d         factor       2           6
e         factor       2           7
f         factor       2           2
g         logical     NA           6
h         integer     NA          12
i         numeric     NA          43
j         factor       2           0
x         factor       2           0
y         factor       2           0

The primary hypothesis is that x can predict y. y is the response variable. I initially removed any observations in which x or y were missing. 192 observations remained for 76 patients. There are missing values, however, for several of the variables: especially h and i, which based on some preliminary tests look like they might be important. There are 50 records in which some variable has a missing value, leaving 142 records with no missing values in any variable.

As far as I understand, glmer() has no way of dealing with missing values and by default will omit the entire record if it has any missing values in the variables being used. I think this may be problematic due to the relatively small sample and effect sizes so I'm not sure what the best way of dealing with this would be.

"id" is the patient id is what I am trying to use as the random effect with all other variables being fixed effects. I've found putting too many variables in the formula will cause glmer to fail to converge. Below are the commands I used for a few simpler models and results:

command:

lr1   <- glmer(formula = y ~ x + 1|id, data = s.dat, family = binomial);

result:

Generalized linear mixed model fit by the Laplace approximation 
Formula: y ~ x + 1 | id 
   Data: s.dat 
   AIC BIC logLik deviance
 238.9 252 -115.4    230.9
Random effects:
 Groups Name        Variance Std.Dev. Corr  
 id     (Intercept) 0.97379  0.98681        
        xpositive   0.01454  0.12058  1.000 
Number of obs: 197, groups: id, 76

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.1343     0.2121  -5.349 8.85e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

command:

lr2   <- glmer(formula = y ~ x * inst + 1|id, data = s.dat, family = binomial);

result:

Generalized linear mixed model fit by the Laplace approximation 
Formula: y ~ x * inst + 1 | id 
   Data: s.dat 
   AIC   BIC logLik deviance
 252.7 288.8 -115.4    230.7
Random effects:
 Groups Name             Variance Std.Dev. Corr                 
 id     (Intercept)      1.016604 1.00827                       
        xpositive        0.072397 0.26907   1.000               
        instI2           1.626269 1.27525  -0.683 -0.683        
        xpositive:instI2 0.135144 0.36762  -0.817 -0.817  0.137 
Number of obs: 197, groups: id, 76

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.1224     0.2099  -5.347 8.96e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

command:

lr3   <- glmer(formula = y ~ x * inst + h + i + 1|id, data = s.dat, family = binomial);

result:

Warning message:
In mer_finalize(ans) : iteration limit reached without convergence (9)

command:

lr4   <- glmer(formula = y ~ x + h + i + 1|id, data = s.dat, family = binomial);

result:

Generalized linear mixed model fit by the Laplace approximation 
Formula: y ~ x + h + i + 1 | id 
   Data: s.dat 
   AIC   BIC logLik deviance
 184.6 217.1 -81.31    162.6
Random effects:
 Groups Name        Variance   Std.Dev.   Corr                 
 id     (Intercept) 1.1663e+01 3.41508554                      
        xpositive   1.2627e+00 1.12369065 -0.284               
        h           7.9415e-02 0.28180626 -0.998  0.229        
        i           2.1064e-07 0.00045896  0.107  0.923 -0.164 
Number of obs: 142, groups: id, 57

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.3359     0.2792  -4.785 1.71e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I don't really understand what this is telling me, or why all my variables are showing up as random effects instead of fixed effects.

Any enlightenment will be greatly appreciated!

Thanks!

Improve this question
$\endgroup$
1
  • $\begingroup$ See my edit below, @ric $\endgroup$
    – Macro
    Jun 19, 2012 at 21:58

1 Answer 1

Reset to default
13
$\begingroup$

Let $Y_i$ be $1$ if patient $i$ is positive and $0$ otherwise and let $X_{i1}, ..., X_{ip}$ be the $p$ predictor variables for patient $i$. A standard tool for determining whether variables are related to an individual's probability that $Y_i = 1$ is the logistic regression model:

$$ \log \left( \frac{ P(Y_i = 1 | X_i) }{P(Y_i = 0| X_i)} \right) = \beta_0 + \beta_1 X_{i1} + ... + \beta_p X_{ip} $$

That is, we fit a regression model to the log-odds that $Y_i = 1$, conditional on the predictors. The coefficients indicate whether or not increases in the predictor increase or decrease the probability that $Y_i = 1$ and, in particular, the exponentiated coefficients can be interpreted as odds ratios. You can then frame your research questions as tests of hypotheses about the regression coefficients. The predictors here can be continuous or categorical, although if they are categorical with more than two categories, you will need to set up dummy variables to look at the effect of the different levels.

You can account for the repeated measurements by doing one of two things, which each have slightly different interpretations:

  • using including random effects that model the correlations between repeated measurements on individuals - in this case the regression coefficients are interpreted as the change in an individual's log-odds that $Y_i = 1$ for a one unit change in the predictor. I'd recommend this approach if the correlations at the individual level are of interest to you. This model can be fit using the lme4 package in R.

  • using GEE which provides robust standard errors for your regression coefficients. This provides unbiased inference despite the fact that there is correlation in the data. The coefficients here are interpreted as the change in the population average log-odds that $Y_i$ for a one unit change in the predictor. I'd recommend this approach if the correlations at the individual level are a nuisance (i.e. not of interest) and you're looking to make inference about the population. This model can be fit using the gee package in R.

  • For more discussion on the difference between random effects models and GEE, see this answer

Edit in response to OP's edit: Good question! This is a subtle issue when you're first learning about it. Based on your comment

"we would want to be able to say when testing new patients in the future that given they are positive for test X, they are Z times more likely to become positive for outcome Y"

it seems that you want to be able to make statements about the change in an individual's odds of a positive test, rather than the population average effect. Therefore, it appears you are more interested in effects at the individual level, not the average effect in the population, so I'd recommend fitting the random effects model.

Regarding your second question, it still could be appropriate to use the logistic model but, from what you've said, centre is an important categorical predictor that should potentially interact with your predictor(s).

Improve this answer
$\endgroup$
9
  • 1
    $\begingroup$ "The coefficients [...] can be interpreted as odds ratios." Do you mean the exponentiated coefficients in the equation $\frac{P(Y_i = 1 | X_i)}{P(Y_i = 0| X_i)} = e^{\beta_0} \cdot e^{\beta_1 X_{i1}} \cdot ... \cdot e^{\beta_p X_{ip}}$? $\endgroup$
    – caracal
    Jun 19, 2012 at 15:38
  • $\begingroup$ Yes, I meant the exponentiated coefficients can be interpreted as odds ratios. Thanks, @caracal. I'll fix it. $\endgroup$
    – Macro
    Jun 19, 2012 at 15:42
  • $\begingroup$ Thanks for the reply! I was under the impression that logistic regression models required continuous predictors so this changes things. $\endgroup$
    – ric
    Jun 19, 2012 at 18:20
  • $\begingroup$ @ric, You're welcome. Like ordinary regression, logistic regression can take both categorical and continuous predictors. Glad I could help. $\endgroup$
    – Macro
    Jun 19, 2012 at 18:20
  • 1
    $\begingroup$ @ric, you're missing some parentheses in your model formula. Try replacing 1|id with (1|id). $\endgroup$
    – Macro
    Jun 22, 2012 at 23:32

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.