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. ~
Assessing linearity
8 messages · Michael Lawrence, Jarrod Hadfield, Ken Knoblauch +2 more
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>:
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. ~
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
Mike Lawrence <Mike.Lawrence at ...> writes:
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
)
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?
Mike
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University
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
Ken Knoblauch Inserm U846 Stem-cell and Brain Research Institute Department of Integrative Neurosciences 18 avenue du Doyen L?pine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.sbri.fr/members/kenneth-knoblauch.html
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, 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>:
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. ~
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
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
Andrew Kosydar, PhD drewdogy at uw.edu On Sun, Oct 24, 2010 at 12:43 PM, Mike Lawrence <Mike.Lawrence at dal.ca> wrote: > 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, >> >> 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>: >> >>> 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. ~ >>> >>> _______________________________________________ >>> R-sig-mixed-models at r-project.org mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models >>> >>> >> >> >> >> -- >> The University of Edinburgh is a charitable body, registered in >> Scotland, with registration number SC005336. >> >> _______________________________________________ >> R-sig-mixed-models at r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models >> > > _______________________________________________ > R-sig-mixed-models at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models >
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:
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, 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>:
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. ~
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
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, 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:
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, 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>:
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. ~
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
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:
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, 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:
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, 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>:
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. ~
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models