mixed model on proportion data
Hello again. This is a query which follows up on the very useful answer provided by Drew Tyre who suggested creating new variables and adding the random term you'll see in the model below. Quick reminder: the analysis aims at understanding how the presence of "bodyguard" ants regulates leaf damage caused by lepidopteran herbivores on a focal plant species. For this, I set up fourteen pairs of plants in a nature reserve and then assigned each plant in each pair to one of two treatments, either "with ants" or "without ants", by applying a physical barrier at the base of the plant. In the following nine weeks I measured four response variables once a week on a subsample of ten randomly chosen leaves of each plant. I thus have nine repeated measures on each plant.
From analyses done over the weekend I know excluding ants from the plants
(variable "treat" in the dataset) has a positive effect on the number of
butterfly eggs and larvae found on plants.
Now I want to know if this results in a reduced damage to the leafs
(variable "dam" in the data set), which is measured as the percentage of
the area eaten by herbivores.
My response variables are:
ants: number of ants (this was measured just to check the physical barrier
had worked OK)
eggs: number of butterfly eggs
larve: number of butterfly larvae
dam: percent leaf damage (percentage eaten by larvae)
My explanatory variables are:
treat: treatment (two levels: "con" and "sin" mean with and without ants,
respectively)
pair: plant pair
date
#I provide a workable example below
require(lme4); require(tidyverse); require(lubridate)
#read data from Google drive
id <- "0Bzd8I1jr8z_iU0h6R0hxaElseDA" # google file ID
gaston <- read.table(sprintf("
https://docs.google.com/uc?id=%s&export=download", id), head=T)
#create variables following suggestion by Drew Tyre
gaston <- mutate(gaston,
plant = paste0(pair,treat),
date = dmy(date),
day = as.numeric((date - min(date))))
#After a few hours of scavenging through books, articles, blogs and R-lists
I ended up with the following model which turns the % damage into a
proportion and then logit transforms it so I can use lmer.
M1 <- lmer(logit(dam/100) ~ treat + day + (1|pair/plant), data=gaston)
summary(M1)
#Am I headed in the right direction? Any advice would be greatly
appreciated.
Best,
Mariano
*Dr. Mariano Devoto*
Profesor Adjunto - C?tedra de Bot?nica General, Facultad de Agronom?a de la
UBA
Investigador Adjunto del CONICET
Av. San Mart?n 4453 - C1417DSE - C. A. de Buenos Aires - Argentina
+5411 4524-8069
http://www.agro.uba.ar/users/mdevoto/