# Strange result in Poisson regression - R

I have performed a Poisson Regression in R, but I got strange results that I cannot find an answer for. My data is like this:

 Aspect_16    Nr_Pereti
1   E         49
2   ENE       73
3   ESE       29
4   N         84
5   NE        77
6   NNE       99
7   NNW       77
8   NW        92
9   S         19
10  SE        20
11  SSE       9
12  SSW       17
13  SW        23
14  W         39
15  WNW       56
16  WSW       25


The Nr_Pereti variable are counts for each level in the 'Aspect_16' column. The model formula and results are:

summary(model_nr_exp)

Call:
glm(formula = tab_gen_exp$Nr_Pereti ~ tab_gen_exp$Aspect_16,
family = poisson)

Deviance Residuals:
[1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                3.8918     0.1429  27.243  < 2e-16 ***
tab_gen_exp$Aspect_16ENE 0.3986 0.1847 2.159 0.030886 * tab_gen_exp$Aspect_16ESE  -0.5245     0.2343  -2.239 0.025169 *
tab_gen_exp$Aspect_16N 0.5390 0.1798 2.998 0.002714 ** tab_gen_exp$Aspect_16NE    0.4520     0.1827   2.473 0.013386 *
tab_gen_exp$Aspect_16NNE 0.7033 0.1747 4.026 5.66e-05 *** tab_gen_exp$Aspect_16NNW   0.4520     0.1827   2.473 0.013386 *
tab_gen_exp$Aspect_16NW 0.6300 0.1769 3.562 0.000368 *** tab_gen_exp$Aspect_16S    -0.9474     0.2703  -3.505 0.000456 ***
tab_gen_exp$Aspect_16SE -0.8961 0.2653 -3.377 0.000733 *** tab_gen_exp$Aspect_16SSE  -1.6946     0.3627  -4.673 2.97e-06 ***
tab_gen_exp$Aspect_16SSW -1.0586 0.2815 -3.761 0.000169 *** tab_gen_exp$Aspect_16SW   -0.7563     0.2528  -2.992 0.002769 **
tab_gen_exp$Aspect_16W -0.2283 0.2146 -1.064 0.287468 tab_gen_exp$Aspect_16WNW   0.1335     0.1956   0.683 0.494845
tab_gen_exp\$Aspect_16WSW  -0.6729     0.2458  -2.738 0.006182 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 2.9506e+02  on 15  degrees of freedom
Residual deviance: 2.1316e-14  on  0  degrees of freedom
AIC: 120.29

Number of Fisher Scoring iterations: 3


I don't have any Deviance Residuals for the model, and when I try to plot the model, it gives me this error:

Error in qqnorm.default(rs, main = main, ylab = ylab23, ylim = ylim, ...) :
y is empty or has only NAs
1: not plotting observations with leverage one:
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning -Inf


What did I do wrong here? Thanks.

You are estimating a fully saturated model, so your model exactly reproduces the counts in your data. So that is why your residuals are all 0.

If you want to plot your model, you can simply plot your data, as the two are exactly the same in this case.

• So does this mean that the model is designed wrong? I am new to the statistics field. – Litwos Jul 19 '15 at 15:53
• No, it just means that your model does not simplify your data. You could make the argument that the purpose of a model is to simplify reality. In that case your "model" is not a model but just another way of printing your dataset. That could be considered wrong or right depending on what you want to do. – Maarten Buis Jul 19 '15 at 15:59
• I am trying to find the influence of the Aspect_16 levels on the Nr_Pereti variable. – Litwos Jul 19 '15 at 16:10

Your problem seems to me to cry out for treatment of your response in terms of the circular or periodic variable aspect using trigonometric functions. Using each aspect as a separate predictor not only produces the problem you met; it is unduly pessimistic.

I did this in Stata; the equivalent should be trivial in R. I know that you can still use glm in R, but you need to convert the compass directions to sine and cosine of aspect.

clear

input obs str3 Aspect_16 Nr_Pereti
1   E         49
2   ENE       73
3   ESE       29
4   N         84
5   NE        77
6   NNE       99
7   NNW       77
8   NW        92
9   S         19
10  SE        20
11  SSE       9
12  SSW       17
13  SW        23
14  W         39
15  WNW       56
16  WSW       25
end

label def asp 1 "N" 2 "NNE" 3 "NE" 4 "ENE" 5 "E" 6 "ESE" 7 "SE" 8 "SSE" 9 "S"
label def asp 10 "SSW" 11 "SW" 12 "WSW" 13 "W" 14 "WNW" 15 "NW" 16 "NNW", add
encode Aspect_16, gen(Aspect) label(asp)
replace Aspect = 22.5 * (Aspect - 1)

sort Aspect
l, sep(4)

gen sine = sin(_pi * Aspect/180)
gen cosine = cos(_pi * Aspect/180)
poisson Nr sine cosine
regplot Aspect, xla(0 "N" 90 "E" 180 "S" 270 "W" 360 "N")


The first half of the code converts the compass direction labels to a numeric equivalent in degrees from North. Then the simplest model to try uses one sine and one cosine term. And it seems to work quite well. Note that a maximum and minimum half a circle apart is an inevitable consequence of one sine and cosine term. Fortunately that also seems natural for many ecological and environmental phenomena.

Any Stata users reading this can download regplot after search regplot.

There is a puff for this kind of modeling here.

• Thanks for the answer. I computed the sin and cos for the Degrees. In the regression model, should the interaction be sin*cos or sin+cos? – Litwos Jul 19 '15 at 16:58
• Interactions are multiplicative in general. In this case, the standard approach is to add terms of the form sin(2 Aspect), cos(2 Aspect), etc. – Nick Cox Jul 19 '15 at 17:04