draft of new mixed models book

I just received this email through the r-sig-mixed-models group, from Douglas Bates:

I have finally gotten the chapter drafts of the book, "lme4:
Mixed-effects Modeling with R", that I am writing to the point where I
think they can be helpful to others.  I believe I have properly
committed them to the R-forge site but they will not appear for an
hour or so.  They should be available under
http://lme4.R-forge.R-project.org/book/

Chapters 1, 2, 4 and 5 are currently available as are the .R files for
those chapters.  Chapters 4 and 5 will have material added to them.  I
don't expect to change chapters 1 and 2 until I get reviewers reports.

When pointing out errors, inconsistencies, bad writing, obtuse
explanations, etc., please be gentle.

lme - lmer mixed models

Dear Kim,

first of all, please use lmer() from library "lme4", instead of lme. I think that is currently better maintained and developed than the "nlme" library.
The model you suggest is a random regression per individual (obs_nr)?
The syntax would become

library(lme4)
cocon_mixed<-lmer(MiWi~host_spec*ins_spec+(date|obs_nr),data=cocon_nash)

It should work provided there are several observations of MiWi for each level of obs_nr, which seems to be the case given your explanation of the data.

Cheers, Tom

Mixed effects models

Hi,

I am currently trying the exercises of Crawley on mixed effects model (ex. 7). I remade the example on page 33/34 and am currently trying tho do this analysis with my own data. I have a response variable (size=MiWi), measured at different dates (date) on the same individual (obs_nr = random variable). Further more I have insects species (ins_spec) and host species (host_spec) as fixed variables. When I run it I get the error message:

cocon_mixed<-lme(fixed=MiWi~host_spec*ins_spec,data=cocon_nash,random=~date|obs_nr) Error in lme.formula(fixed = MiWi ~ host_spec * ins_spec, data = cocon_nash, :
nlminb problem, convergence error code = 1
message = iteration limit reached without convergence (9)

What is going wrong?

Kim

Course GLM Groningen

Hi,
I just returned from the glm course in Groningen. I will invite all participants to this blog in the coming days. Meanwhile, I already add a link to a page with several simple scripts:

http://tomvandooren.eu/tuto.html

As you will see, these files haven't been updated at all.
If things in them have become obsolete, please let me know.

Cheers, Tom

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

multinom

Hi Tessa,
To compare the fit of two models, you can first fit each one separately, and compare them using the anova() function. It will do a likelihood ratio test.
I think you have to construct your response differently, and not use "replica" to group observations.
The response variable has to be a matrix, with a row for each replica, and the counts on each patch in the columns. The explanatory variables will have as many entries as there are lines in the matrix. I now see that your covariates are patch-specific, which is more involved than what multinom() allows. Multinom() only allows an offset with that structure.
Maybe you will have to fit your data using a log-linear poisson model then, or function vglm() from the "VGAM" library can also handle multinomial data. Agresti, Categorical Data Analysis has a lot of info on the analysis of multinomial data.
Scripts for that book are available at http://www-stat.stanford.edu/~owen/courses/306a/Splusdiscrete2.pdf.
Cheers, Tom

planned comparisons

Hi Nicola,
you told me it's data on selection experiments. Are the data rally three blocks, meaning crosses between selection lines started from three different base populations?

If you write it like this:
model<-aov(egg~block+block:cross)
there there are no main effects for cross, and you only want to compare crosses within blocks. That means that you only want contrasts between different groups in the same block.

To compare these groups, first look at the coefficients of the model, and what they represent.
Reformulating the model like this probably makes comparisons easier:
model<-aov(egg~block:cross)
because all estimates are means per level.

With four crosses per block, are the coefficients of block 1, the first four coefficients in the model?
Say, you want to make simultaneous confidence intervals for all pairwise comparisons for this block:

# contrast matrix

cmatrix<-matrix(rep(0,6*4),nrow=6)
cmatrix[1,1]<-1;cmatrix[1,2]<--1;
cmatrix[2,1]<-1;cmatrix[2,3]<--1;
cmatrix[3,1]<-1;cmatrix[3,4]<--1;
cmatrix[4,2]<-1;cmatrix[4,3]<--1;
cmatrix[5,2]<-1;cmatrix[5,4]<--1;
cmatrix[6,3]<-1;cmatrix[6,4]<--1;

cmatrix

# when this matrix is multiplied with the vector of coefficients,
# every element of the result, is an estimated difference between two means.

# library(multcomp)

# psubset selects the parameters of the block you want to investigate.

simint(model,psubset=1:4,cmatrix=cmatrix)

You should probably make a contrast matrix for all 18 pairwise comparisons simultaneously.

cmatrix<-matrix(rep(0,3*6*4*3),nrow=3*6)
cmatrix[1,1]<-1;cmatrix[1,2]<--1;
cmatrix[2,1]<-1;cmatrix[2,3]<--1;
cmatrix[3,1]<-1;cmatrix[3,4]<--1;
cmatrix[4,2]<-1;cmatrix[4,3]<--1;
cmatrix[5,2]<-1;cmatrix[5,4]<--1;
cmatrix[6,3]<-1;cmatrix[6,4]<--1;

cmatrix[7,5]<-1;cmatrix[7,6]<--1;
cmatrix[8,5]<-1;cmatrix[8,7]<--1;
cmatrix[9,5]<-1;cmatrix[9,8]<--1;
cmatrix[10,6]<-1;cmatrix[10,7]<--1;
cmatrix[11,6]<-1;cmatrix[11,8]<--1;
cmatrix[12,7]<-1;cmatrix[12,8]<--1;

cmatrix[13,9]<-1;cmatrix[13,10]<--1;
cmatrix[14,9]<-1;cmatrix[14,11]<--1;
cmatrix[15,9]<-1;cmatrix[15,12]<--1;
cmatrix[16,10]<-1;cmatrix[16,11]<--1;
cmatrix[17,10]<-1;cmatrix[17,12]<--1;
cmatrix[18,11]<-1;cmatrix[18,12]<--1;

simint(model,psubset=1:12,cmatrix=cmatrix)

Multinomial

Hi Tessa,


Answers:

1) how do I define my distribution in a glm when I have multinomial data (5 patches)?

Well, the glm() function cannot handle multinomial data, but function multinom() from the "nnet" library can.
Another package which fits a multinomial probit model is "MNP", there is also the "multinomrob" library. I think multinom() is easiest to use.


2) the same question but for a generalized linear mixed model (I would like to include time as a factor)?

When you need to add a random effect, "MNP" is the only option in R as far as I can see. Things become complicated in that case. Research on multinomial random effects models is still ongoing.
Some references: Multinomial logit random effects models, Statistical Modelling (2001) (J. Hartzel, A. Agresti, and B. Caffo); The analysis of ordered categorical data: An overview and a survey of recent developments, invited discussion paper for the Spanish Statistical Journal, TEST (2005) (I. Liu and A. Agresti).

3) In Crawley they talk about GAM’s (general additive models). What I understand is that for these models you don’t need to define a distribution because they use non parametric smoothers. Can these models be used for multinomial data?

Generalized additive models use non-parametric smoothers to model a non-linear dependence of a response variable on a predictor variable.
The distribution of the response variable is among the list of usual suspects: binomial, poisson, normal, gamma.
It has been claimed that package "gbm" (R) or "MART" (S-Plus) can do this for multinomial data, but I haven't tried those yet.

4) Spatial distributions are difficult. Package "spatstat" could contain the solution. I am not very familiar with analysing spatial patterns.

random effects

Hi Fransesc,
you seem to have only a small number of levels for the plot effects.
The consequence will be that estimates of variance between these random effects are very imprecise, and that it is not very solid to state that these levels are a representative sample from a large population of plot levels.
I would fit all effects as fixed.
Continuous variables such as time cannot define levels random effects, but they can appear in random regressions. Anyway, if you want to fit higher-order models of time, all you have to do is make new variables, like
timesq<-time^2
and add them to the model.
In a model with interactions of these time variables with treatment (or other factors), different polynomials of time will be fit for each treatment group.
Hope that helps,
Tom

variable lengths

Hi Cas,
something has gone wrong with attaching the data, or calculating new variables from old ones. The error message you get means that different variables in your model have different lenghts (so there's a mismatch).
for example, for a model of the form
glm(y~x+z)
you can check the lengths of x,z and y by
length(y)
length(x)
length(z)
The numbers returned should all be the same for a correct model specification. That makes it easy to spot where things went wrong.
good luck, Tom

mixed model

Hi Barbara,
I guess you want to fit a linear mixed model, nlme is specific for non-linear models.

Did you try this code, which assumes a normally distributed response:

model1<-lme(mfreq~avergTemp+altitude, random=~1|country)
summary(model1)

or this code for poisson response:

library(MASS)
model2<-glmmPQL(mfreq~avergTemp+altitude, random=~1|country, family=poisson)
summary(model2)

library(lme4)
lmer(mfreq~avergTemp+altitude+(1|country), family=poisson, method="Laplace")

but first you should check whether country is a factor!
summary(country)

Cheers, Tom

count data and dispersion

First of all,
long ago I picked up the idea that corrections for overdispersion are only necessary when the scale parameters is above 3 to 4. The idea is understandable, you should only correct if likely to be present, but the exact reason for the "3 to 4" limit is now unclear to me.
I tried to locate where I got the idea from, and it seems to be from here. I did not find any confirmation of that rule of thumb anywhere else.
However, Lindsey (1999) suggests, based on the analysis of examples, that corrections are necessary when the overdispersion parameter is at least two.

You should keep in mind, however, that overdispersion is impossible in some cases (so you should then never correct for it):

- when the dependent variable is a Bernoulli 0 - 1 variable
- when the maximal model considered, is the saturated model

Overdispersion can also be caused by

- using the wrong link function
- a missing covariate
- the necessity to transform covariates
- outliers

so you should check these potential causes.

Concerning underdispersion, Venables and Ripley (2002) show that dispersion estimates indicating underdispersion can be caused by small counts. See e.g. MASS4 p.208.
The ratio of residual deviance to degrees of freedom can be seriously biased downward, for extreme p (and small n) in the case of binomial data, and for low lambda.

Instead of relying on rules of thumb, one can also model distributions for counts which allow both under-and overdispersion.
Some possibilities are the double binomial (poisson) distribution and the multiplicative binomial (poisson) distribution. Using AIC, you can then incorporate the decision to correct for under(over)dispersion into your model selection procedure.
You can do glm with these distributions in R, using the "gnlm" library by (Jim Lindsey. It is not easy to work with.
I have adapted some examples from Lindsey 1999 Models for Repeated Measurements to indicate how you can fit such models.
The R input code is here, and the data file here.

GLM course

The glm course has started this morning. Thomas Tully, Wolf Mooij and Tom Van Dooren are the instructors. All participants have been invited to become members of this blog.
Tom

Quickstart R

I've made a one-hour tutorial to get you started in R statistical software. It's stored at:
http://biology.leidenuniv.nl/~dooren/Rglimpse.html
Enjoy!
Tom