[R-sig-ME] Resulat from GLMM, error bars for predictions

[Quentin Schorpp]

Hello Everyone,

I'm using glmmadmb() to model abundance data (counts) of soil organisms 
(e.g. earthworms).  My design covers agricultural fields of five 
different age classes . Every age-class has three field replicates. 
Additionally every field was sampled 3times during the investigation 
(atumn, spring, autumn -> sampling campaign).
with 4 samples taken randomly at each field (1 Sample = 0.25 m²). 
Several environmental parameters were assessed for each field but never 
for one of the four samples explicitly.
Hence the environmental data is often redundant for the samples, 
especially when climatic measurements were even similar for more ththan 
one field.

My hypothesis is that abundance increases with the age_class

My model is:

model <- glmmadmb(abundance ~ age_class + samplingCampaign + 
environmental1 + env2 + env3 + (1|field), data, family="poisson")

age_class = ordered factor
sampling campaign = continous, difference to the first sampling in days 
(first sampling always 0)
env(1-n) = continous
total number of samples = 180

Dispersion factor is 1.45

I do model validation with
1. plot pearson residuals against fitted values
2. plot pearson residuals against each covariate in the model
3. make a histogram of the residuals

In my opinion everything looks ok.

Now I have the really really big problem: *I just don't know how to 
present the results.*
I'd like to do a barplot with mean abundances for age_class and standard 
errors and do Post Hoc tukey test to look at differences between the 
factor levels. But i just don't know how to to these Post-Hoc tests.

I've got one approach for extracting predictions and standard errors for 
predictions using a test dataset with mean environmental variables:

test.data=expand.grid(age_class=levels(data$age_class),
                                         samplingCampaign = data$samcam),
                                         env1 = mean(data$env1),
                                         env2 = max(data$env2))

pred.abundance <- cbind(test.data,
                                               predict(model, test.data, 
type="link", se.fit=TRUE),
                                               abundance.response = 
predict(model, test.data, type="response"))

pred.anc <- within(pred.abundance, {
                                          anc <- 4*exp(fit)
                                          LL <- 4*exp(fit - 1.96 * se.fit)
                                          UL <- 4*exp(fit + 1.96 * 
se.fit)  })

Then I make a plot and get INCREDIBLY large standard errors and in 
contrast to the boxplot of the predicitons 
(plot(data$age_class,predict(model, type="response")), the abundance is 
not increasing with the age_class. I multiply by 4 since i want to 
present the results per m²

Do you know where the mistake is?

I would appreciate if you could help me with this analysis, since I'm 
trying to learn GLMM for more than a year and i can't ask a real person 
here at this Institution. Thanks in advance,
Quentin


-- 

Quentin Schorpp, M.Sc.
Thünen Institute of Biodiversity
Bundesallee 50
38116 Braunschweig (Germany)

Tel:  +49 531 596-2524
Fax:  +49 531 596-2599
Mail: quentin.schorpp at ti.bund.de
Web:  http://www.ti.bund.de

The Johann Heinrich von Thünen Institute, Federal Research Institute for Rural Areas, Forestry and Fisheries – Thünen Institute in brief –
consists of 15 specialized institutes that carry out research and provide policy advice in the fields of economy, ecology and technology.


	[[alternative HTML version deleted]]