We sample plots of high elevation rocky outcrops, counting the number of American Pika within each plot; plot sizes are the same size.
Poisson Regression
We model the counts of American pika (\(y_{i}\)) at each site \(i = 1...n\) as a Poisson random variable with parameter \(\lambda\) being a function of our \(p\) site-level variables in \(n\) x \(p\) design matrix (\(\textbf{X}\)) and coefficients \(\boldsymbol{\beta}\) as,
model =glm(y~site,data=dat, family =poisson(link ='log') )
summary(model)
Call:
glm(formula = y ~ site, family = poisson(link = "log"), data = dat)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6549 0.1140 5.747 9.09e-09 ***
siteSite 2 -1.0116 0.2207 -4.584 4.56e-06 ***
siteSite 3 0.1221 0.1565 0.780 0.435
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 341.36 on 119 degrees of freedom
Residual deviance: 305.76 on 117 degrees of freedom
AIC: 496.33
Number of Fisher Scoring iterations: 6
Numerical Optimization of Likelihood
model$converged
[1] TRUE
model$method
[1] "glm.fit"
model$boundary
[1] FALSE
model$iter
[1] 6
model$control
$epsilon
[1] 1e-08
$maxit
[1] 25
$trace
[1] FALSE
Side-Bar
Negative Log-Likelihood Function
#Here is our negative log-likelihood function with three#parameters - beta0, beta1, and beta2 (1)#inputs = design matrix X neg.log.like =function(par,X) {#linear model of parameters par and design matrix (X) lam=par[1]*X[,1]+par[2]*X[,2]+par[3]*X[,3]#neg log-likelihoodsum(-dpois(y,lambda =exp(lam),log =TRUE)) }
Numerical Optim. for MLEs
#Use optim function with initial values and define the lower and upper limits of the possible values fit1 <-optim(par =c(0,0,0), #startX = X,fn = neg.log.like,method ="L-BFGS-B",lower =c(-10, -10, -10),upper =c(10, 10, 10) )
X =model.matrix(~Independent Variables)#all variables in dat (additive)X =model.matrix(~., data = dat)#all variables in dat (pair-wise interactions)X =model.matrix(~.^2, data = dat)#all variables in dat (three way- interactions)X =model.matrix(~.^3, data = dat) model =glm(y~0+ X, ...)
Poisson Regression (offset)
What if our counts of pika are at plots with different sizes?
Plot size (\(\textbf{A}\)) needs to be controlled for. But, we don’t want to estimate an effect as it’s part of the design. Rather, we want to model the rate - counts per unit area.
Additive & Interaction combinations of categorical and continuous variables
site (categorical)
herb.cover (continuous)
dist.trail (continuous)
head(dat)
y site herb.cover dist.trail site.re area
1 0 Site 1 -0.22630093 0.06679655 Site 1 2
2 1 Site 2 0.02221148 -0.24641488 Site 2 5
3 0 Site 3 -1.01517375 -0.73612867 Site 3 2
4 1 Site 1 1.18038309 0.07586011 Site 1 3
5 0 Site 2 -2.04470642 -0.81595409 Site 2 2
6 2 Site 3 2.44433686 0.17929858 Site 3 2
