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?
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 pearson 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.
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.
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!