Significantly more sightings over the years?

This is probably something very easy for many of you, but I am not sure if I am doing this right. I am a biologist and have been looking at the occurrence (per year) of a certain species in a certain area from 2001 to 2016. As this is a rare species, in two years we saw none, in the other years we saw 1-10 individuals. It seems like the animals are being seen more in recent years and I think there is a significant increase of them over the years. But how do I test this? I have simplified my dataset to just years and a number of sightings in that year (of different individuals, so no double counts): If I run a linear regression in SPSS my p-value is significant. But is this the right way to do it? • Yes, sure! Is it better like this? – Rockhouse May 16 '17 at 14:12
• It looks like time series analysis would be better., – Michael Chernick May 16 '17 at 14:19
• Rare count can be modelled by Poisson distribution, so a poisson regression might do ... or time series, but with a few sightings of a rare animal I doubt if time correlation is important. More sightings might be because of more effort, are you sure the effort was approximately constant? (effort might be used as an ofset in a poisson model, so if you have data about effort, post that also) – kjetil b halvorsen May 16 '17 at 14:27
• @kjetilbhalvorsen Unfortunately I do not have effort data, but I do think effort was fairly constant. Do you think this dataset is too small for a statistical conclusion? – Rockhouse May 16 '17 at 14:30

As kjetil writes, there is rather little data here - too little to draw truly firm conclusions.

My personal impulse in such a situation is to look at multiple models and see whether they agree in principle.

Let's first draw a graph.

year <- 2001:2016
sightings <- ts(c(1,0,3,2,0,4,2,2,5,4,4,10,8,3,6,5),start=2001)
plot(sightings,type="o",pch=19) Well, the graph does look rather convincing, even if it is only 16 data points. I'd certainly rather bet on the 2017 observation to be larger than 4, rather than 4 or less.

Here is the R analogue to the linear regression you calculated in SPSS. Note that the p values match:

summary(lm(sightings~year))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -796.7588   220.6300  -3.611  0.00283 **
year           0.3985     0.1098   3.628  0.00274 **

As kjetil suggested, a Poisson regression would also make sense:

summary(glm(sightings~year,family="poisson"))

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -227.69850   61.49336  -3.703 0.000213 ***
year           0.11395    0.03058   3.726 0.000194 ***

This model also finds a significant trend. However, the linear and the Poisson regression posit very specific trends: a linear one in the first case and an exponential one in the second case.

An alternative would be to test the correlation between sightings and year. The correlation again calculates a linear trend:

cor.test(sightings,year,method="pearson")

t = 3.628, df = 14, p-value = 0.002742

Note how the p value is exactly the same as for the simple linear model.

Alternatives would be the and the correlations, which test for any kind of monotonous relationship between years and sightings, not only linear or exponential ones:

cor.test(sightings,year,method="kendall")

z = 2.9576, p-value = 0.0031

cor.test(sightings,year,method="spearman")

S = 156.75, p-value = 0.0004915

Although both tests output warnings, because they cannot calculate exact p values if ties are present, they again find significant trends.

Finally, as Michael Chernick notes, you actually have a time series, so a time series analysis might be useful. Your count data really call for an INAR model or similar, but there are really no common count data time series models that account for trend, so I'll just fit an ARIMA model and an ETS one:

library(forecast)
auto.arima(sightings)

Series: sightings
ARIMA(0,1,0)

sigma^2 estimated as 8:  log likelihood=-36.88
AIC=75.76   AICc=76.07   BIC=76.47

ets(sightings)

ETS(A,N,N)

Call:
ets(y = sightings)

Smoothing parameters:
alpha = 0.3891

Initial states:
l = 1.3

sigma:  2.3277

We note that auto.arima() models an ARIMA(0,1,0) process, which means that it believes that first differences are white noise. First differences again indicate a trend. Finally, ETS is the only one that does not find a trend, it only finds additive error (the first "A"), no trend ("N") and no seasonality ("N"). However, it finds a very large smoothing value of $\alpha = 0.39$, so it thinks your sightings might be a weak kind of a random walk. Note that these models are fitted using information criteria, so it doesn't make sense to assign a p value to the trends they find (or not).

In summary, most of these different models do find a trend, even if they are designed to look at different kinds of trends (linear, exponential, general monotone, first differences). This, together with the plot, would certainly be enough to convince me that there is indeed a trend in your data.

• Thank you so much for your effort and explanation. That really helped a lot! I will have to look at all of this again in detail but it makes a lot more sense now. – Rockhouse May 17 '17 at 8:03

Small counts like with a few sightings of a rare animal are usually well modelled by a Poisson distribution, so Poisson regression is indicated. Variability in number of sightings of a rare animal could well be caused by a varying level of effort, below I assume this is not the case, that is, effort is approximately constant over time. (If you have effort data, it could be included as an offset via the argument offset=log(effort) in the glm command below).

I show below an analysis of your data using Poisson regression in R, which I think is more appropriate than the analysis you presents in the question. R code is:

year      <-  2001:2016
sighting  <-  c(1,0,3,2,0,4,2,2,5,4,4,10,8,3,6,5)
dat  <-  data.frame(year,sighting)
# Poisson regression model:
mod.0  <-  glm(sighting ~ year, family=poisson, data=dat)
summary(mod.0)
confint(mod.0)
plot(mod.0)

Which produces the following output:

> mod.0  <-  glm(sighting ~ year, family=poisson, data=dat)
> summary(mod.0)

Call:
glm(formula = sighting ~ year, family = poisson, data = dat)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.0802  -0.7293  -0.2194   0.8220   2.0665

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -227.69850   61.49336  -3.703 0.000213 ***
year           0.11395    0.03058   3.726 0.000194 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 33.801  on 15  degrees of freedom
Residual deviance: 18.770  on 14  degrees of freedom
AIC: 67.106

Number of Fisher Scoring iterations: 5

> confint(mod.0)
Waiting for profiling to be done...
2.5 %       97.5 %
(Intercept) -351.90183933 -109.9909922
year           0.05539646    0.1757022

(I do not show here the plots, which do not indicate any problems with the model)

Note that the confidence interval on the slope parameter is far away from zero, so some increase in abundance with time is indicated. With so few data effective model criticism (based on residuals and such) is not really possible, so you need to believe in the model!

• +1. The final discussion is especially important, since these data exhibit an interesting characteristic. Because sightings are assigned to years, rather than specific dates, we might ask how sensitive the results are to slight changes in the dates. The answer is that the most recent three years are noticeably more influential on the coefficient estimates than any of the previous years. Thus, if the purpose is to estimate rates of change over time, we should note that (1) this Poisson model posits exponential increase and (2) the last three years of data begin to challenge that assumption. – whuber May 16 '17 at 15:39
• Thank you very much for the effort and explanation! I have to look at this again but it does make more sense than my SPSS try. – Rockhouse May 17 '17 at 8:05