random effects

Dear Tom,
my apologies for the long delay between you request and this answer.
First of all, why not include precipitation ("Precip") as a main effect too?

Using libray "nlme" and lme():
Random intercepts per station:
random= ~1|station
Don't forget the vertical bar!
Including random slope:
random= ~Precip|station
or
random= ~1+Precip|station

Make sure that all stations have different names, otherwise you run the risk that stations with identical names at different distances will be considered as a single station.

You could also consider to use the updated library now called "lme4".
Then the random effects are simply included as extra terms in the model statement.

Cheers, Tom

random components of lme-model

Hi everyone,

I'm currently analysing a large dataset with measurements of nutrients in the Dutch marine systems in R. I have a variable "OffshoreDist" representing the distance from the sampling location to the coast. I also calculated the total amount of precipitation over The Netherlands and called this variable "Precip".What I want to do is link mean annual dissolved inorganic nitrogen (DIN) to the amount of precipitation. I expect that this relationship is dependent on the offshore distance. Distance classes can contain 1 or more sampling stations.

That's why I propose these two models with different random components:
model1 <- lme(log(DIN) ~ Precip:factor(OffshoreDist) + factor(OffshoreDist),
random= ~1station)
model2 <- lme(log(DIN) ~ Precip:factor(OffshoreDist) + factor(OffshoreDist),
random= ~Precipstation)

Are these the right models for random intercepts per station (model1) and random slope+intercept per station (model2)? Or should this be something like:

random = ~ factor(OffshoreDist) -1 station (model 1)
and
random = ~ Precip:factor(OffshoreDist)station (model 2)

I've read Pinheiro & Bates, 2000 which is a good reference, but explanation of the formulae only by means of examples does not seem to work for me.

Thanks in advance.

Cheers,

Tom Van Engeland

moved to new blogger

Hi,
I've moved the blog into the new google template. Everyone logging in will be invited to move her/his account too. A nice aspect of the change is that the blog can now be searched easily!
Cheers, Tom

F lmer

Hi Marie,
do you really need to use the quasibinomial distribution? Is there evidence of overdispersion in the maximal model you fit?

If not, then you can use the binomial distribution.
That allows you to use likelihood-ratio tests for the fixed and random effects.
using
anova(model,modelsimplified,test="Chisq")

If there is clear overdispersion, then I guess the best thing to do would be resampling according to the model, if that works for the quasibinomial "distribution". You need to use mcmcsamp() for that.
Otherwise, it will become hand-work to get some of those hotly debated F-ratios out....
Cheers, Tom

Longitudinal Data

Finally, the participants in the course Longitudinal and Incomplete Data have been invited. Transcripts for R of the SAS code presented during the course will be added to the weblog soon.
Cheers, Tom

Loglinear models

Hi Nicola,
the best manual is http://www.statslab.cam.ac.uk/~pat/Splusdiscrete2.pdf, based on the book by Agresti.
You can basically use glm with a poisson link for log-linear models, I guess it is also discussed in the Venables and Ripley book.
Good luck, Tom

lmer, random effect with four levels

Dear Amir,
the "place" random effect in your dataset has only four levels, and they are nested within plants. The error which is returned states that the variance covariance matrix is singular, meaning that variation in intercept, slope and their covariance cannot be fitted. If you fit the random effect as (1|place), you will read that the random effect variance is effectively zero.
It does not make sense to include place both as fixed and as random, as in you second to fourth models. Since the number of levels of the random effects is very low, you could even just model them as fixed, using glm. The idea of random effects is that they are a representative sample from a large population, which is hard to defend with just two levels sampled per group!
More importantly, the distribution of your data does not look poisson at al. There are many zero observations. I would convert the counts to presence/absence data and analyse those as bernoulli probabilities. That works to some extent:

b<-ifelse(Predators.Allspecies.>0,1,0)
model<-lmer(b~julian.date+offset(log.time)+(1|place),family=binomial)

I've added you models to your post, so other visitors to this blog can also read and respond to your question,
all the best, Tom