This week I am writing about how to estimate count data models in R. Similarly to last week, I am using data from Sohngen et al. (2000). The authors asked a sample of recreationists for the number of past single-day trips to Lake Eirie beaches. The data I am using can be downloaded here, under Lake Eirie Beach Data.
Data can come in many forms. While in methods such as contingent valuation data are discrete, in travel cost method applications we usually analyze data that is visitation frequency. This data is usually in the format of a count, ranging from zero (meaning no visits) to the maximum number of trips of an avid recreationist.
To take a look at the Lake Eirie visitation data, I can use the summary function to see the descriptive statistics for the number of trips variable.
> summary(lake_erie_beach_data$tthsbch2) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000 3.000 5.000 6.542 11.000 15.000
The mean number of trips is around 6.5 per respondent. However, the median is slightly lower, which means that the distribution of trips is skewed to the left.
In order to look at the distribution of trips, I will plot a histogram:
Indeed, there are a lot of respondents who have been to Lake Eirie beaches less than 4 times.
In environmental economics, count data models are commonly used when applying the travel cost method. Whereas a model estimated by OLS estimates the marginal effect of the explanatory variables on the number of trips, count data models estimate how the probability of a given count changes with changes to those explanatory variables.
The advantage of count data models is that they do not predict negative values for the count variable. In this previous post, I obtained some negative values for the predicted number of trips when using an OLS regression. This does not happen when using count data models instead.
The Poisson model
There are two commonly used count data models: the Poisson and Negative Binomial models. However, there are other options when it comes to count data models. Which model should be used?
The Poisson model is perhaps the most commonly used count data model. When we use the Poisson regression to explain trip frequency, we assume that the number of trips has a Poisson distribution. We can then model linearly the logarithm of the expected number of trips given some relevant explanatory variables.
Here are the explanatory variables in our example:
- trcst30, i.e. the total travel cost to get to the Lake Eirie beach.
- Male, which takes the value 1 if the repondent is male, and 0 if the respondent is female.
- Age, which represents the age of the respondent.
- People, i.e. number of peole in the group.
- Income, i.e. the household’s annual income.
- Destin, which takes the value 1 if the beach was the trip’s primary destination, or 0 otherwise.
- Pbeach, i.e. the proportion of time spent on the beach.
- Zcctcst3, zcltcst3 and zgvtcst3, which are total travel costs to three alternative state sites.
In R, we can run a Poisson regression using the glm function. The terminology is very similar to a simple linear regression (with the function lm): the dependent variable tthsbch2 (number of trips) is before the ~ sign, followed by all the explanatory variables. Also, I identify the dataframe that has all these variables: data=lake_erie_beach_data. The difference is that, when calling the glm function, I need to specify the family=poisson.
model_p <- glm(tthsbch2 ~ trcst30 + male + Age + people + income +
destin + pbeach + zcctcst3 +
zcltcst3 + zgvtcst3,
family="poisson", data=lake_erie_beach_data)
And the output looks like the following:
> summary(model_p) Call: glm(formula = tthsbch2 ~ male + Age + people + income + destin + pbeach + zcctcst3 + zcltcst3 + zgvtcst3, family = "poisson", data = lake_erie_beach_data) Deviance Residuals: Min 1Q Median 3Q Max -3.8786 -1.7063 -0.3689 1.2110 6.3101 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.643e+00 3.145e-01 5.224 1.75e-07 *** trcst30 -1.531e-02 1.950e-03 -7.850 4.16e-15 *** male 2.345e-01 4.605e-02 5.092 3.55e-07 *** Age 3.858e-03 1.593e-03 2.422 0.0154 * people -3.752e-02 9.283e-03 -4.042 5.31e-05 *** income 3.867e-06 3.590e-06 1.077 0.2815 destin 4.551e-01 6.774e-02 6.719 1.83e-11 *** pbeach -5.021e-02 6.990e-02 -0.718 0.4725 zcctcst3 -5.881e-04 4.672e-03 -0.126 0.8998 zcltcst3 3.955e-03 4.068e-03 0.972 0.3308 zgvtcst3 -5.580e-03 3.000e-03 -1.860 0.0629 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 1327.1 on 309 degrees of freedom Residual deviance: 1004.2 on 300 degrees of freedom (35 observations deleted due to missingness) AIC: 2066.8 Number of Fisher Scoring iterations: 6
We can conclude by looking at the results that the number of trips is expected to decrease with increases in travel cost and number of people in the group, for example. The number of beach trips is expected to increase with marginal increases in the age of the respondent and whether getting to the beach was the primary destination of the trip.
The Negative Binomial Model
One of the underlying assumptions of the Poisson model is that the conditional mean of the dependent variable equals its conditional variance.
The unconditional means and variance of the Lake Eirie is quite different. The unconditional mean is:
mean(lake_erie_beach_data$tthsbch2) [1] 6.515942
The standard deviation is:
> sd(lake_erie_beach_data$tthsbch2) [1] 5.22621
The unconditional variance is two times the standard deviation, so it is around 25.
It turns out that many times this assumption fails. That is, the conditional variance exceeds the conditional mean. We can check if this also holds true in our data by typing:
pr <- residuals(model_p,"pearson") phi <- sum(pr^2)/df.residual(model_p) phi
phi has a value of 4. This means that the conditional variance is four times higher than the conditional mean. Hence, it would be more appropriate to use other count data models, rather than the Poisson
This also happens in other studies. Blaine et al. (2015) reported that the average number of trips per visitor was around 21 (unconditional mean), whereas the standard deviation was around 30. In another study (Roussel et al., 2016), the average number of visits per year was 7.15, whereas the standard deviation was much higher (19.44).
If this happens, i.e. conditional variance exceeds the mean, then we say the data is over-dispersed. In such a case, the standard errors of the Poisson model and likely wrong and we can use the Negative Binomial model instead.
The Negative Binomial model is a generalization of the Poisson model. In R, it can be estimated using the package MASS. Within MASS, I can use the glm.nb function, whose terminology is as simple as using the lm function for linear regressions.
library(MASS)
model_nb <- glm.nb(tthsbch2 ~ trcst30 + male + Age + people + income +
destin + pbeach + zcctcst3 +
zcltcst3 + zgvtcst3,
data=lake_erie_beach_data)
The estimated model is the following:
Estimate Std. Error z value Pr(>|z|) (Intercept) 1.239e+00 5.010e-01 2.473 0.013395 * trcst30 -1.435e-02 3.109e-03 -4.615 3.94e-06 *** male 2.217e-01 8.554e-02 2.592 0.009549 ** Age 2.527e-03 3.040e-03 0.831 0.405770 people -1.727e-02 8.292e-03 -2.082 0.037314 * income -4.683e-09 5.651e-06 -0.001 0.999339 destin 4.236e-01 1.193e-01 3.550 0.000385 *** pbeach -2.098e-04 1.370e-01 -0.002 0.998778 zcctcst3 4.841e-03 7.359e-03 0.658 0.510671 zcltcst3 1.152e-03 6.907e-03 0.167 0.867581 zgvtcst3 -2.954e-03 4.792e-03 -0.616 0.537640
The estimated coefficients are very similar to those estimated with the Poisson model. However, many of the Negative Binomial coefficients have lost statistical significance (e.g. Age). This usually happens, i.e. the standard errors estimated with the Negative Binomial model are broader that those estimated with the Poisson model.
References:
Blaine, T. W., Lichtkoppler, F. R., Bader, T. J., Hartman, T. J., & Lucente, J. E. (2015). An examination of sources of sensitivity of consumer surplus estimates in travel cost models. Journal of environmental management, 151, 427-436.
Roussel, S., Salles, J. M., & Tardieu, L. (2016). Recreation demand analysis of sensitive natural areas from an on-site survey. Revue dEconomie Regionale Urbaine, (2), 355-384.
2 Comments