Skip to content

Assessing linearity

8 messages · Michael Lawrence, Jarrod Hadfield, Ken Knoblauch +2 more

#
Hi folks,

I have developmental data collected across several grades (1-6). I
would like to be able to assess whether there are any linear or
non-linear trends across grade. Does it make sense to run a first lmer
treating grade as continuous, obtain the residuals, then run a second
lmer treating grade as a factor? That is:

fit1 = lmer(
    formula = response ~ (1|individual)+grade_as_numeric
    , data = my_data
    , family = gaussian
)
my_data$resid = residuals(fit1)
fit2 = lmer(
    formula = resid ~ (1|individual)+grade_as_factor
    , data = my_data
    , family = gaussian
)


As I understand it, fit1 will tell me if there are any linear trends
in the data, while fit2 will tell me if there are any non-linear
trends in the data in addition to the linear trends obtained in fit1.

If this is sensible, how might I apply it to a second binomial
response variable given that the residuals from a binomial model are
not 0/1?

Cheers,

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tr.im/mikes_public_calendar

~ Certainty is folly... I think. ~
#
Hi Mike,

You would be better off trying out something like polynomials or  
splines, For example:

  fit1 = lmer(
      formula = response ~ (1|individual)+poly(grade_as_numeric,n),
      , data = my_data
      , family = gaussian
  )

where n is the order of the polynomial. n=1 would fit the same model  
as your original fit1, although the covariate (and the regression  
parameter) would be scaled by some number. When n=6 the model would be  
a reparameterised version of your model fit2. When 1<n<6 you would be  
working with a non-linear relationship in between these two extremes,  
although the model is still linear in the parameters.

Cheers,

Jarrod








Quoting Mike Lawrence <Mike.Lawrence at dal.ca>:

  
    
#
Mike Lawrence <Mike.Lawrence at ...> writes:
Or just make grade and ordered factor and see whether there
are significant higher order terms.

If the response variable is binomial, why not use glmer with a
binomial family?

Ken
#
Ah, this looks promising! So how does this sound:

I typically assess the evidence for a relationship between the
predictor and response variables by comparing the AIC values for a
model including the predictor to a model without it. In the case of
grade_as_numeric, I'd do:

fit_null = lmer(
    formula = response ~ (1|individual)
    data = my_data
)
fit_null_AIC = AIC(fit_null)

fit_alt = lmer(
    formula = response ~ (1|individual) + grade_as_numeric
    data = my_data
)
fit_alt_AIC = AIC(fit_alt)

grade_loglikratio = fit_null_AIC - fit_alt_AIC

Now, if I wanted to check whether there is a quadratic component to
the grade effect, I'd first compute an analogous likelihood ratio for
the quadratic fit compared to the null:
fit_alt_quad = lmer(
    formula = response ~ (1|individual) + poly(grade_as_numeric)^2
    data = my_data
)
fit_alt_quad_AIC = AIC(fit_alt_quad)
grade_quad_loglikratio = fit_null_AIC - fit_alt_quad_AIC

Then compute a final log likelihood ratio between the improvement over
the null caused by grade versus the improvement over the null caused
by grade as a quadratic:

grade_lin_vs_quad_LLR = grade_loglikratio - grade_quad_loglikratio

I could repeat this for higher-order polynomials of grade, each
compared to the order directly below it, to develop a series of
likelihoood ratios that describe the relative improvement of the fit.

Does this sound appropriate?

Cheers,

Mike

--
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tr.im/mikes_public_calendar

~ Certainty is folly... I think. ~
On Sun, Oct 24, 2010 at 6:41 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
#
Hi Mike,

Sounds as though you are debating about model comparison techniques.
I would suggest using multimodel inference if you do not have any one
model that is strongly supported by AIC values.  Burnham & Anderson
published a book on the technique and they have two papers (2000 &
2002) that lay out the theoretical framework and methodology pretty
clearly.

Good luck!

Andrew
#
Hi Mike,

Couple of things.

poly() is used like this

a <- rnorm(100)
b <- rnorm(100)

m2 <- lm(a~poly(b,2))
summary(m2)

To compare a model with a linear predictor and a quadratic polynomial
you can do this:
m0 <- lm(a~1)
m1 <- lm(a~poly(b,1))
m2 <- lm(a~poly(b,2))
anova(m1,m2)

In the case of likelihood fitted models you'll get a likelihood ratio test.


You wouldn't normally fit a quadratic term instead of the linear term,
rather in addition, so m2 has 1 extra parameter. You can compare
multiple models at once with anova(m0,m1,m2), or you can go the AIC
route and calculate deltaAIC for your set of candidate models.

Unless you expect the non-linearity to be polynomial shaped you might
be better off fitting splines with varying number of knots/df.



andydolman at gmail.com
On 24 October 2010 18:43, Mike Lawrence <Mike.Lawrence at dal.ca> wrote:
#
Thanks Andrew. Jarrod mentioned using splines as well, and I do think
they sound promising since I indeed do not have strong a priori
reasons to expect polynomial shaped trends.

I'm uncertain, however, how to choose values for the "degree" and "df"
arguments to the bs. I see that df has to be equal-to-or greater than
degree, and less than the number of unique values in grade minus 2
(minus one makes it equivalent to factorizing grade, which I'm doing
already). I also see that degree=1 and df=1 is equivalent to the
linear case, which I'm already estimating in the original lmer fit.

However, that still leaves a large parameter space (df at degree:
2 at 1;3 at 1,4 at 1, 2 at 2, 3 at 2; 4 at 2; 3 at 3; 4 at 3; 4 at 4), and I'm not sure it's
appropriate to compare the linear case to fits from every point in
this space. Thoughts?

Mike
On Mon, Oct 25, 2010 at 12:59 PM, Andrew Dolman <andydolman at gmail.com> wrote:
#
Hi Mike,

I don't see why you would want to fit such a complicated trend to your
data. If you need anything more than a cubic then just treat grade as
a factor. I would start by fitting a model with natural splines ns(),
not bs() with 1 and 2 degrees of freedom and then see whether your 2
degree model fits better - this would be a reasonable test for
non-linearity in grade. If you expect the effect of grade to jump
around all over the place then treating it as a factor if probably
more appropriate. Even a cubic is probably overkill, do you expect
performance to ever decline with grade?



andydolman at gmail.com
On 25 October 2010 20:15, Mike Lawrence <Mike.Lawrence at dal.ca> wrote: