Dear list members,

I would appreciate some advice on how to model nonlinearity in my dataset.

The dataset contains >1000 observations on variables Q and M. The
observations (actually species) are partitioned into around 50 groups. Each
group spans only a small part of the total range of M. The principle
research focus is on the relationship between Q and M, and whether this
relationship displays heterogenity among groups. This is relatively easy to
explore using:

m1 <- lmer( Q ~ M +  (1|g) )
m2 <- lmer( Q ~ M +  (M|g) )

So far so good. However, I also want to test some recent theories that Q~M
should vary with M, particularly that the relationship between slope and M
might be negative curvilinear.

Until recently, I thought that it would be appropriate to extract the
conditional modes from the model and seek correlations separately. However,
recent posts by Reinhold Kliegl on this list have made it clear that this
would not be appropriate.

So I've been wondered how to model this in lmer. I was thinking of fitting
2nd order polynomials to model the curvilinearity:

m3 <- lmer( Q ~ M + I((M+20)^2) +  (M|g) ) #M is centered on it's mean, so I
need to add 20 to make all values positive
m4 <- lmer( Q ~ M + I((M+20)^2) +  (M + I((M+20)^2)|g) ) )
m5 <- lmer( Q ~ M + I((M+20)^2) + (M|g) + (I((M+20)^2)|g) ) )

I'm not sure how I would interpret differences in the fit of these models,
and I'm sure that other (probably better) solutions exist.

It has also been proposed that variance in the Q~M relationship should
decline with increasing M. Is this a question about heteroscedasticity that
I need to model explicitly (the residuals of m2 are quite well-behaved)?

I would be extremely grateful for any advice on this topic.

Best wishes, Nick

	[[alternative HTML version deleted]]


