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?
 A: 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.)
A: 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
