I am think that it is possible to analyse a model with just random effects but I am not sure as I have never done it. I am looking for guidance on whether it is appropriate, what assumptions I need to be aware of, and how to do it properly.

From my study of an insect;

  • I have a response variable (age at death, "age")
  • Two treatments ("Treat1" and "Treat2") both of which have two levels (Treat1 has "A" and "B", and Treat2 has "P" and "Q")
  • There is also 40 genotypes (1-40)
  • With four replicates (w,x,y,z) of each combination of Genotype/Treat1/Treat2
  • Each replicate contains 50 individuals

Put simply, my data looks like 32000 rows of this:

Treat1  Treat2  Genotype  Block  Individual   Age   
A       P       1         w      1            23
A       P       1         w      2            35
A       P       1         w      3            44
.       .       .         .      .            .
.       .       .         .      .            .
.       .       .         .      .            .
B       Q       40        z      50           76     

I would like to know if each combination of Treat1 and Treat2 (AP,AQ,BP,BQ) have genetic genetic variation - i.e. is there variation between my 40 genotypes within each treatment combination?

I think I need a model for each of AP, AQ, BP, and BQ, along the lines of

Age ~ Genotype [ Treat1 == "A" & Treat2 == "P"] * Block [ Treat1 == "A" & Treat2 == "P"]

Where Genotype and Block are random effects. I hear Gamma distribtions are better to use in lifespan (time to death) models.

My questions are:

a. Is this an appropriate way to show whether or not my genotypes have variation?

b. Can I build the four models as defined above or is that a really poor way of doing it?

c. If possible, what functions should I be using in R (lm, glm, lmer... & summary, summary.lm, aov, anova...)?

d. What should I expect, if gamma is more suitable than gaussian, to see when I compare plot(model) for gamma compared to gaussian?

This is currently my model...

AP= df$Treat1=="A" & df$Treat2=="P"
apmodel<- lmer(df$Age[AP]~(1|df$Genotype[AP])+(1|df$Block[AP]))

Which I think is right but I'm not sure what to do with the output..

> summary(apmodel)
Linear mixed model fit by REML 
Formula: df$Age[AP] ~ (1 | df$Genotype[AP]) + (1 | df$Block[AP]) 
       AIC   BIC logLik deviance REMLdev
     57343 57371 -28667    57336   57335
    Random effects:
     Groups           Name        Variance Std.Dev.
     df$Genotype[AP]  (Intercept) 17.23798 4.15186 
     df$Block[AP]     (Intercept)  0.15416 0.39263 
     Residual                     93.18777 9.65338 
    Number of obs: 7757, groups: df$line[AP], 40; df$Block[AP], 4

Fixed effects:
            Estimate Std. Error t value
(Intercept)  49.9948     0.6939   72.05

Is there genetic variance??


There are several things you could do to test whether there is genetic variance.

First, however, I wonder why you want separate models for "each combination of Treat1 and Treat2 (AP,AQ,BP,BQ)"? I don't know anything about the substantive area of application here, and I may be misunderstanding your data, but I think you can have 5 varying intercepts / random effects here: Block, Genotype, Treat1, Treat2, and Treat1.Treat2 -- an interaction term that you create and add to the dataset. One model for everything.

Anyway, back to your question.

First, a formal test of significance can be conducted by running this model using ML (using REML = FALSE):

apmodel1 <- lmer(Age ~ (1|Genotype) + (1|Block), df, REML = FALSE)

Then running a model without the Genotype effect:

apmodel2 <- lmer(Age ~ (1|Block), df, REML = FALSE)

And performing a likelihood ratio test:

anova(apmodel2, apmodel)

Second, and more informally, you might calculate statistics such as the intraclass correlation coefficient (ICC) to provide a measure of how much variance Genotype is accounting for. Try the ICC.lme function from the psychometric library.

Finally, the most interesting method would be to produce plots of predicted effects, with their corresponding uncertainty estimates for each Genotype, Treatment, etc. The focus here shifts from Genotype as a factor, to each level of Genotype as a (modelled / varying) effect.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.