# Appropriate method for determining difference between means?

I'm looking at the calling rate of bird species when prey availability is altered. So there are two groups, a control (no manipulation) and experimental group (where prey has been increased). I'm interested in looking to see if there is a difference in mean calling rate between two different groups but would also like to include other factors (like environmental factors) and see how they might influence the calling rate.

I am using R for my analysis and am wondering if anyone has any insight on what types of statistical tests I might want to use?

• How big are the groups? Nov 5, 2014 at 23:40
• If you have counts (of calls) for a given duration, I would model the rate of calling by a Poisson GLM with an offset (see ?glm; the offset would be logged duration). You can include treatments and covariates of interest (environmental variables, prey availability, etc.) as predictors. Temporal and site effects (for example) can be included as random effects (see ?glmer in package lme4). Nov 6, 2014 at 2:21
• Thanks for the response! I have 200 surveys at the control stations and 173 surveys at the experimental surveys. The surveys were 20 minutes long so its the number of calls heard within that 20 minute period. I didn't measure the actual prey available but the experimental survey stations had more prey available than the control. Nov 6, 2014 at 19:52
• @gung OK, I've turned the comment into an answer. Nov 10, 2014 at 1:01
• @NatePope, nice job. Nov 10, 2014 at 3:33

I'd model the total number of calls over a survey duration as a Poisson GLM with an offset.

What is an offset? Let $y_i$ be the observed number of calls and $\lambda_i$ be the expected number of calls, and $t_i$ be the duration (in minutes) of the $i$th survey. Let $X$ be a matrix of regressors and $\beta$ be the associated vector of regression coefficients.

A model which does not include duration would look like this:

\begin{aligned} y_i \sim \mathrm{Poisson}(\lambda_i) \\ \log(\lambda_i) = x_i^T\beta \end{aligned}

The above equation expresses the total calling rate as a function of covariates. However, if the survey duration differs among samples, this model will be incorrect. Instead, we'd want:

\begin{aligned} y_i \sim \mathrm{Poisson}(\lambda_i/t_i) \\ \log(\lambda_i/t_i) = x_i^T\beta \end{aligned}

Which would instead model the per-minute calling rate as a function of covariates (essentially standardizing the expected number of calls by the amount of sampling effort). Rearranging:

\begin{aligned} \log(\lambda_i/t_i) = \log(\lambda_i) - \log(t_i) = x_i^T\beta \\ \Rightarrow \log(\lambda_i) = x_i^T\beta + \log(t_i) \end{aligned}

$\log(t_i)$ is referred to as an 'offset' or as 'exposure'. The offset won't really matter for inference if every survey has the same duration, but I'll consider an example where this isn't the case, to provide a more general answer.

Using your sample size, let's say there is an average calling rate of 6 calls per minute in the control and 4 calls per minute in the experimental units. Further, let's pretend that surveys are split among seven sites.

As an example environmental covariate, we'll also suppose that vegetation cover influences calling rate such that a 1% increase in vegetation cover causes a multiplicative increase of 0.001 in calling rate per minute.

Some R code to provide an example:

## simulated data
set.seed(101)
treatment   <- rep(c(0, 1), c(200, 173)) # dummy variable for treatment
n.site      <- 7 # number of sites
site.effect <- rnorm(7, 0, 0.4) # site deviations
site.labels <- sample(1:7, length(treatment), replace = T) # random site labels
#   for surveys
intercept   <- log(6) # rate per minute in control
exp.effect  <- log(6)-log(4) # diff. between control and experimental unit,
#   in calling rate per minute
veg.cover   <- rbeta(length(treatment), 0.3*7, 0.7*7)*100 # fake vegetation cover data
veg.effect  <- 0.001 # multiplicative effect of vegetation on calling rate
l.duration  <- log(rpois(length(treatment),20)) # logged survey duration (random)
lambda      <- exp(intercept + site.effect[site.labels] + treatment*exp.effect +
veg.cover*veg.effect + l.duration) # linear predictor
y           <- rpois(length(treatment), lambda) # simulated response variable
#   (ie. total number of calls)
## fit model to data
# first ignoring site effects,
m0 <- glm(y ~ treatment + veg.cover + offset(l.duration), family = poisson)
# now with random site effects
library(lme4)
m1 <- glmer(y ~ treatment + veg.cover + offset(l.duration) + (1|site.labels),
family = poisson) # ignore warnings

summary(m1)
## we recover accurate estimates for treatment effect, vegetation cover,
## and a less accurate estimate of the variance of the site effect.
## note that the regression coefficients are on a log (ie. a multiplicative) scale


The function offset() forces the regression coefficient to 1.

You may also want to consider overdispersion (I typically include overdispersion as a matter of course). One way to do so is to include an individual-level random effect (essentially, a multiplicative residual) in the call to glmer():

id <- 1:length(treatment)
glmer(y ~ treatment + veg.cover + offset(l.duration) + (1|site.labels) + (1|id),
family = poisson)


The variance for the individual-level random effect is small in this example because we do not include overdispersion in the simulated response. However, it is very very often the case that actual count data will be overdispersed--in other words, there will be additional variation in the response which is attributable to unobserved factors (ie. in your case, unmeasured factors may include additional environmental variables, biotic factors, population density, observer effects, etc.)

• Great answer. Can glm of base R be used for this? What would be the command to use glm to analyze this?
– rnso
Nov 10, 2014 at 1:16
• @rnso offset() works within glm(). I gave an example of this in the code snip. Overdispersion with glm() can be modeled either by using a quasi- distribution (for example, family=quasipoisson) or with glm.nb() from MASS. Nov 10, 2014 at 1:40
• Thanks for the thorough response! Greatly appreciate it. My surveys were all 20 minutes long so the duration is the same for all the data. So I assume I wouldn't need to include the offset in this case? Nov 10, 2014 at 19:40
• @pprest00 That's correct. You can verify that the inference will not change, by fitting a model with a constant offset, and comparing against a model without that offset. In the example I give above, only the parameter estimate for the intercept would change. Nov 12, 2014 at 0:35

Assuming $N$ is reasonably large (maybe 35 or more per group), then you can use a $t$-test, which assumes the variance is both groups is equally large.

This makes sense because you randomized, so we would not expect the variance of one the observations to be systematically different across groups.

To run this in R, see:

http://www.stat.columbia.edu/~martin/W2024/R2.pdf

• Thanks for the response. I ran a t-test to see if there was a difference in calling but that doesn't take into account the other environmental factors that may be confounding the results. So i was interested in how to include these as well. Nov 6, 2014 at 19:53
• Then the commenter is right about using GLM with Poisson. I thought you were A/B testing. Nov 6, 2014 at 23:09