Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Introduction:
When we deal with data that has a response variable of integer type, using a linear regression may violate the normality assumption and hence all the classical statistic tests would fail to evaluate the model. However, as we do with logistic regression models, the generalized linear model GLM can be used instead here by specifying the suitable distribution.
The possible distributions for this type of data are the discrete distributions poisson and negative binomial. The former is the best choice if the mean and the variance of the response variable are closer to each other, if they are not however and we persist using this distribution we may cause the rise of the overdispersion problem of the residuals. As a solution thus, we can use the latter distribution that does not have this restriction.
There is another alternative if neither the poisson distribution nor the negative binomial are suitable called the Quasi maximum likelihood. The advantage of this method is that uses only the relationship between the mean and the variance and does not require any prespecified distribution. Moreover, its estimators are approximately as efficient as the maximum likelihood estimators.
Data preparation
To well understand how to model the count data we are going be using Doctorvisits data from AER package, in which the variable visits will be our target variable, so let’s call this data with the packages that we need along this article.
ssh <- suppressPackageStartupMessages
ssh(library(performance))
ssh(library(ModelMetrics))
ssh(library(corrr))
ssh(library(purrr))
ssh(library(MASS))
ssh(library(tidyverse))
ssh(library(AER))
ssh(library(broom))
data("DoctorVisits")
doc <- DoctorVisits
glimpse(doc)
Observations: 5,190
Variables: 12
$ visits <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, ...
$ gender <fct> female, female, male, male, male, female, female, female,...
$ age <dbl> 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.1...
$ income <dbl> 0.55, 0.45, 0.90, 0.15, 0.45, 0.35, 0.55, 0.15, 0.65, 0.1...
$ illness <dbl> 1, 1, 3, 1, 2, 5, 4, 3, 2, 1, 1, 2, 3, 4, 3, 2, 1, 1, 1, ...
$ reduced <dbl> 4, 2, 0, 0, 5, 1, 0, 0, 0, 0, 0, 0, 13, 7, 1, 0, 0, 1, 0,...
$ health <dbl> 1, 1, 0, 0, 1, 9, 2, 6, 5, 0, 0, 2, 1, 6, 0, 7, 5, 0, 0, ...
$ private <fct> yes, yes, no, no, no, no, no, no, yes, yes, no, no, no, n...
$ freepoor <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, n...
$ freerepat <fct> no, no, no, no, no, no, no, no, no, no, no, yes, no, no, ...
$ nchronic <fct> no, no, no, no, yes, yes, no, no, no, no, no, no, yes, ye...
$ lchronic <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, n...
This data from Australian health survey where visits is the number of doctor visits in past two weeks with 11 features listed above.
First we list the summary of the data to inspect any unwanted issue.
summary(doc)
visits gender age income
Min. :0.0000 male :2488 Min. :0.1900 Min. :0.0000
1st Qu.:0.0000 female:2702 1st Qu.:0.2200 1st Qu.:0.2500
Median :0.0000 Median :0.3200 Median :0.5500
Mean :0.3017 Mean :0.4064 Mean :0.5832
3rd Qu.:0.0000 3rd Qu.:0.6200 3rd Qu.:0.9000
Max. :9.0000 Max. :0.7200 Max. :1.5000
illness reduced health private freepoor
Min. :0.000 Min. : 0.0000 Min. : 0.000 no :2892 no :4968
1st Qu.:0.000 1st Qu.: 0.0000 1st Qu.: 0.000 yes:2298 yes: 222
Median :1.000 Median : 0.0000 Median : 0.000
Mean :1.432 Mean : 0.8619 Mean : 1.218
3rd Qu.:2.000 3rd Qu.: 0.0000 3rd Qu.: 2.000
Max. :5.000 Max. :14.0000 Max. :12.000
freerepat nchronic lchronic
no :4099 no :3098 no :4585
yes:1091 yes:2092 yes: 605
As we see we do not have missing values and the visits values ranges from 0 to 9 but it should be of integer type rather than double. Similarly, the variable illness should be converted to factor type since it has a few different values.
doc$visits<-as.integer(doc$visits) doc$illness <- as.factor(doc$illness) tab <- table(doc$visits) tab 0 1 2 3 4 5 6 7 8 9 4141 782 174 30 24 9 12 12 5 1
The best thing we do to start analyzing the data is by displaying the correlation coefficient of each pair variables we have. Thus, any particular predictor that has high correlation with the target variable could be highly likely to be relevant in our future model. Notice that our target variable is not continuous hence we will use the spearman correlation. As required by correlate function from corrr package, all the variables must be of numeric type so we convert all the factor to integer.
doc1 <-modify_if(doc, is.factor, as.integer)
notice that we have stored the result in another object doc1 to keep save our original data.
M <- correlate(doc1, method="spearman")
rplot(shave(M), colours=c("red", "white", "blue" ))+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Looking at this plot all the correlations has low values. however, these correlations assess only the monotonic relations, they say nothing about any other form of relation.
First let’s compare the empirical distribution of the variable visits and the theoretical poisson distribution with \(\lambda\) equals the visits mean mean(doc$visits), and the total number of observations is 5190.
pos <- dpois(0:9,0.302)*5190
both <- numeric(20)
both[1:20 %% 2 != 0] <- tab
both[1:20 %% 2 == 0] <- pos
labels<-character(20)
labels[1:20 %% 2==0]<-as.character(0:9)
barplot(both,col=rep(c("red","yellow"),10),names=labels)
As we see the two distributions are more or less closer to each other. Let’s now check the negative binomial distribution by first estimate the clumping parameter \(k=\frac{\bar x^2}{s^2-\bar x}\).
k<-mean(doc$visits)^2/(var(doc$visits)-mean(doc$visits))
bin<-dnbinom(0:9,0.27,mu=0.302)*5190
both1<-numeric(20)
both1[1:20 %% 2 != 0]<-tab
both1[1:20 %% 2 == 0]<-bin
labels<-character(20)
labels[1:20 %% 2==0]<-as.character(0:9)
barplot(both1,col=rep(c("red","yellow"),10),names=labels)
With this distribution it seems that the empiricall distribution is more closer to the negative binomial than the poisson distribution.
Note: This data has very large number of zeros for the outcome compared to the other values which means that any trained model that does not take into account this anomaly will be biased to predict more likely the zero value. However, at the end of this article I will show two famous models to handel this type of count data called Haurdle model and zero_inflated model.
Data partition
In oreder to evaluate our model we held out 20% of the data as testing set.
set.seed(123) index<-sample(2,nrow(doc),replace = TRUE,p=c(.8,.2)) train<-doc[index==1,] test<-doc[index==2,]
Poisson model
This model belongs to the generalized linear model families, so in the function glm we set the argument family to poisson. In practice this model is sufficient with a wide range of count data.
set.seed(123) model1<-glm(visits~., data=train, family ="poisson") tidy(model1) # A tibble: 16 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -2.70 0.141 -19.2 9.14e- 82 2 genderfemale 0.193 0.0620 3.11 1.88e- 3 3 age 0.436 0.184 2.37 1.77e- 2 4 income -0.161 0.0928 -1.74 8.23e- 2 5 illness1 0.944 0.113 8.35 6.76e- 17 6 illness2 1.21 0.118 10.3 1.15e- 24 7 illness3 1.11 0.132 8.43 3.51e- 17 8 illness4 1.28 0.140 9.13 7.14e- 20 9 illness5 1.44 0.139 10.4 2.34e- 25 10 reduced 0.126 0.00560 22.6 6.85e-113 11 health 0.0348 0.0112 3.10 1.91e- 3 12 privateyes 0.111 0.0795 1.39 1.64e- 1 13 freepooryes -0.344 0.190 -1.81 7.00e- 2 14 freerepatyes 0.0377 0.104 0.363 7.16e- 1 15 nchronicyes 0.0186 0.0732 0.254 7.99e- 1 16 lchronicyes 0.0255 0.0916 0.279 7.81e- 1
As we see all the variables are significant except for the income so we remove this variable and reestimate again.
set.seed(123) model1<-glm(visits~.-income, data=train, family ="poisson") tidy(model1) # A tibble: 15 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -2.83 0.121 -23.4 7.36e-121 2 genderfemale 0.213 0.0609 3.51 4.56e- 4 3 age 0.479 0.183 2.62 8.91e- 3 4 illness1 0.946 0.113 8.38 5.44e- 17 5 illness2 1.21 0.118 10.3 8.29e- 25 6 illness3 1.12 0.132 8.50 1.93e- 17 7 illness4 1.28 0.140 9.17 4.71e- 20 8 illness5 1.45 0.139 10.5 1.05e- 25 9 reduced 0.126 0.00560 22.6 1.12e-112 10 health 0.0350 0.0112 3.11 1.84e- 3 11 privateyes 0.100 0.0793 1.27 2.06e- 1 12 freepooryes -0.290 0.188 -1.55 1.22e- 1 13 freerepatyes 0.0683 0.102 0.667 5.05e- 1 14 nchronicyes 0.0171 0.0731 0.235 8.15e- 1 15 lchronicyes 0.0282 0.0914 0.308 7.58e- 1
For the interpretation of the coefficient estimates, we should exponentiate these values to get the marginal effect since the poisson model uses the log link function to preclude negative values. For continuous predictor, say age, if this predictor increases by one year, ceteris-paribus, we expect the doctor visits will be \(exp(0.47876624)=1.614082\) times larger. whereas, for categorical predictor, say gender, the female has \(exp(0.21342446)=1.23791\) larger doctor visits than male.
By looking at p-values all the predictors are significant. However, we have to check other statistics and metrics.
glance(model1)
# A tibble: 1 x 7
null.deviance df.null logLik AIC BIC deviance df.residual
<dbl> <int> <dbl> <dbl> <dbl> <dbl> <int>
1 4565. 4154 -2685. 5399. 5494. 3486. 4140
since the deviance value 3485.905 is lower than the degrees of freedom 4140, we will then worry about overdispersion problem. Fortunateley, the AER package provides a super easy way to test the significance of this difference via the function dispersiontest.
dispersiontest(model1)
Overdispersion test
data: model1
z = 6.278, p-value = 1.714e-10
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
1.397176
If our target variable really follows poisson distribution then its variance \(V\) should be approximately equal to its mean \(\mu\), which is the null hypothesis of the following dispersiontest test against the alternative hypothesis that the variance of the form: \[V=\mu+\alpha.trafo(\mu)\]
Where the trafo is an hyperparameter that should be specified as an argument of this test. The popular choices for this argument are:
- trafo = NULL (default): \(V=(1+\alpha)\mu\)
- trafo = 1: \(V=\mu+\alpha.\mu\)
- trafo = 2: \(V=\mu+\alpha.\mu^2\)
For the first choice if true, then the data will be better modeled by quasi-poisson model than poisson model.
For the last ones if one of them is true then the negative binomial will be better than poisson model.
Now once the trafo is defined the test estimates \(\alpha\), such that:
- if \(\alpha = 0\) : equidispersion (The null hypothesis)
- if \(\alpha < 0\) : underdispersion
- if \(\alpha > 0\) : overdispersion
Therefore, the result of the test will depend on the direction of the test, where we have two.sided, greater (default) for the overdispersion, and less for underdispersion.
With this in mind the output of the above test (with the default values) tested the overdispersion against the quasi-poisson model, and since the p-value is very small 1.714e-10 then we have overdispersion problem, suggesting the use of quasi-poisson model instead.
Now let’s test the negative binomial now.
dispersiontest(model1, trafo = 1)
Overdispersion test
data: model1
z = 6.278, p-value = 1.714e-10
alternative hypothesis: true alpha is greater than 0
sample estimates:
alpha
0.3971763
The test suggested the use of negative binomial with linear function for the variance with very tiny p-value 1.714e-10. This model is known as NB1 (with linear variance function).
dispersiontest(model1, trafo = 2)
Overdispersion test
data: model1
z = 7.4723, p-value = 3.939e-14
alternative hypothesis: true alpha is greater than 0
sample estimates:
alpha
0.95488
If the relation is in quadratic form then this model is called NB2. And since this p-value 3.939e-14 is smaller than the previous one then NB2 could be more appropriate than NB1.
Quasi poisson model
The first test suggested the use of quasi-poisson model, so let’s train this model with the same predictors as the previous one.
set.seed(123) model2<-glm(visits~.-income, data=train, family ="quasipoisson") tidy(model2) # A tibble: 15 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -2.83 0.140 -20.2 1.25e-86 2 genderfemale 0.213 0.0705 3.03 2.47e- 3 3 age 0.479 0.212 2.26 2.39e- 2 4 illness1 0.946 0.131 7.24 5.36e-13 5 illness2 1.21 0.136 8.89 9.15e-19 6 illness3 1.12 0.152 7.34 2.49e-13 7 illness4 1.28 0.162 7.92 2.91e-15 8 illness5 1.45 0.160 9.06 2.01e-19 9 reduced 0.126 0.00647 19.5 4.79e-81 10 health 0.0350 0.0130 2.69 7.14e- 3 11 privateyes 0.100 0.0917 1.09 2.74e- 1 12 freepooryes -0.290 0.217 -1.34 1.81e- 1 13 freerepatyes 0.0683 0.119 0.576 5.64e- 1 14 nchronicyes 0.0171 0.0846 0.203 8.39e- 1 15 lchronicyes 0.0282 0.106 0.266 7.90e- 1
This model uses the quasi maximum likelihood which gives the same coefficient estimates but with different (corrected) standard errors.
Since here also all the variables are significant We see that the models are the same except the correction of the standard errors which are now more larger. In other words, the poisson distribution under overdispersion underestimates the standard errors and hence the t test would be biased towards the rejection of the null hypothesis.
To better understand what is going on with quasi-poisson model let’s put the estimates and the standard errors of both models into one table, and we add a column that resulted from dividing the second standard errors vector by the first one.
D1 <- tidy(model1)
colnames(D1) <- NULL
D2 <- tidy(model2)
colnames(D2) <- NULL
cbind(term=D1[,1], estimate1=D1[,2], std1=D1[,3],estimate2=D2[,2], std2=D2[,3]) %>%
mutate(dispersion= std2/std1)
term estimate1 std1 estimate2 std2 dispersion
1 (Intercept) -2.82868582 0.121003941 -2.82868582 0.140023371 1.15718
2 genderfemale 0.21342446 0.060883445 0.21342446 0.070453120 1.15718
3 age 0.47876624 0.183054884 0.47876624 0.211827495 1.15718
4 illness1 0.94643852 0.112984298 0.94643852 0.130743198 1.15718
5 illness2 1.21007304 0.117661029 1.21007304 0.136155018 1.15718
6 illness3 1.11941845 0.131726935 1.11941845 0.152431807 1.15718
7 illness4 1.28108306 0.139695671 1.28108306 0.161653071 1.15718
8 illness5 1.45198517 0.138529459 1.45198517 0.160303554 1.15718
9 reduced 0.12621548 0.005595111 0.12621548 0.006474552 1.15718
10 health 0.03497647 0.011228772 0.03497647 0.012993714 1.15718
11 privateyes 0.10033918 0.079266327 0.10033918 0.091725427 1.15718
12 freepooryes -0.29008562 0.187570838 -0.29008562 0.217053268 1.15718
13 freerepatyes 0.06829904 0.102422323 0.06829904 0.118521089 1.15718
14 nchronicyes 0.01714304 0.073082122 0.01714304 0.084569188 1.15718
15 lchronicyes 0.02818992 0.091425522 0.02818992 0.105795808 1.15718
Note: The first two columns are for the model1, and the last one are for the model 2. Not surprisingly that the result of the last column is constant since this is exactely what the quasi maximum likelihood does, it computes the corrected standard errors from the original ones as follows \(std2=dispersion*std1\), with the dispersion value being estimated as 1.15718. if you want to know where this value came from, the answer is simple. this model computes the sigma of the standardized residuals resulted from the original model. we can thus get this value by specifying the argument type to pear then computing sigma by hand as follows:
resid <- resid(model1, type = "pear") sqrt(sum(resid^2)/4140) [1] 1.15718
Now to test the prediction qualities of our models we use the testing set test by ploting the original and the predicted values. Let’s start by the model1
pred<- predict.glm(model1,newdata=test[test$visits!=0,],type = "response") plot(test$visits[test$visits!=0],type = "b",col="red") lines(round(pred),col="blue")
If you noticed, and due to the large number of zero’s of the target variable, i have intentionally removed all theses values in order to get clearer plot. From this plot we can say that the model does not fit well the data especially the larger values that are not well captured, however this may due to the fact that the data are very skewed towards zero.
To compare different models we can use the root mean-square error and mean absolute error (all the data with zero’s included).
Note: Here we are using the rmse function from ModelMetrics that expects the inpute to be two vectors, and not that with the same name from the performance package that expects the input to be a model object . To avoid thus any such ambiguity you should type this command ModelMetrics::rmse.
pred <- predict.glm(model1, newdata = test, type = "response") rmsemodelp <- ModelMetrics::rmse(test$visits,round(pred)) maemodelp <- mae(test$visits,round(pred)) rmsemodelp [1] 0.7381921 maemodelp [1] 0.284058
By the same way, Now let’s evaluate the quasi-poisson model.
predq<- predict.glm(model2,newdata=test[test$visits!=0,],type = "response") plot(test$visits[test$visits!=0],type = "b",col="red") lines(round(predq),col="blue")
This plot does not seem to be very different from the previous plot. The rmse and mae for this model are computed as follows.
predq <- predict.glm(model2,newdata=test, type = "response") rmsemodelqp <- ModelMetrics::rmse(test$visits,round(predq)) maemodelqp <- mae(test$visits,round(predq)) rmsemodelqp [1] 0.7381921 maemodelqp [1] 0.284058
we will not compare this two models until we finish with all the incoming models and we compare all the models at once.
Negative binomial model
The negative binomial distribution is used as an alternative for the poisson distribution under overdispersion problem.
set.seed(123)
model3<-glm.nb(visits~.-income, data=train)
summary(model3)
Call:
glm.nb(formula = visits ~ . - income, data = train, init.theta = 0.9715923611,
link = log)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8413 -0.6894 -0.5335 -0.3540 3.6726
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.894820 0.135760 -21.323 < 2e-16 ***
genderfemale 0.258968 0.075352 3.437 0.000589 ***
age 0.511867 0.230297 2.223 0.026240 *
illness1 0.880644 0.123264 7.144 9.04e-13 ***
illness2 1.171615 0.130240 8.996 < 2e-16 ***
illness3 1.118067 0.149032 7.502 6.28e-14 ***
illness4 1.263367 0.165370 7.640 2.18e-14 ***
illness5 1.378166 0.169907 8.111 5.01e-16 ***
reduced 0.141389 0.008184 17.275 < 2e-16 ***
health 0.041364 0.015029 2.752 0.005918 **
privateyes 0.086188 0.095173 0.906 0.365149
freepooryes -0.375471 0.223857 -1.677 0.093487 .
freerepatyes 0.144928 0.127751 1.134 0.256602
nchronicyes 0.022111 0.087590 0.252 0.800705
lchronicyes 0.091622 0.114965 0.797 0.425477
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(0.9716) family taken to be 1)
Null deviance: 3208.0 on 4154 degrees of freedom
Residual deviance: 2431.2 on 4140 degrees of freedom
AIC: 5159.6
Number of Fisher Scoring iterations: 1
Theta: 0.972
Std. Err.: 0.103
2 x log-likelihood: -5127.587
As before we visualize the performance of this model as follows.
prednb<- predict.glm(model3,newdata=test[test$visits!=0,],type = "response") plot(test$visits[test$visits!=0],type = "b",col="red") lines(round(prednb),col="blue")
Again this plot also seems to be the same as the previous ones, so to figure out which model is best we use statistic metrics.
prednb<- predict.glm(model3,newdata=test,type = "response") rmsemodelnb<-ModelMetrics::rmse(test$visits,round(prednb)) maemodelnb<-mae(test$visits,round(prednb)) rmsemodelnb [1] 0.7808085 maemodelnb [1] 0.2966184
we will use these ouputs further.
Hurdle model
Originally proposed by Mullahy (1986) this model can take into account the fact that the data has more zeros and also can handle the overdispersion problem. It has two components (or steps), truncated count component defined by the chosen discrete distribution such as poisson or negative binomial, and a hurdle components models zero vs larger counts (that uses censored count distribution or binomial model). In other words, this models asumes that two population distributions underlying the data, one distribution for zero values, and another different distribution the psotive values. For more detail about hurdle and zero inflated models click here
To perform this model we make use of the function hurdle from the package pscl.
hurdle model with poisson distribution.
This model works in two steps. In the first step it uses binary classification to discriminate between the zero values and the positive values, and in the second step uses the traditional (poisson or binomial model, and here we use poisson model) model for positive values.
library(pscl)
set.seed(123)
modelhp<-hurdle(visits~. -income, data=train,dist = "poisson")
summary(modelhp)
Call:
hurdle(formula = visits ~ . - income, data = train, dist = "poisson")
Pearson residuals:
Min 1Q Median 3Q Max
-1.5464 -0.4686 -0.3306 -0.2075 11.0887
Count model coefficients (truncated poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.977535 0.261835 -3.733 0.000189 ***
genderfemale 0.073326 0.098034 0.748 0.454480
age 0.032762 0.287991 0.114 0.909427
illness1 0.370071 0.251920 1.469 0.141833
illness2 0.403514 0.256363 1.574 0.115489
illness3 0.201724 0.278757 0.724 0.469277
illness4 0.420285 0.277573 1.514 0.129990
illness5 0.762209 0.269809 2.825 0.004728 **
reduced 0.111640 0.007967 14.013 < 2e-16 ***
health 0.007682 0.016452 0.467 0.640554
privateyes -0.215649 0.129860 -1.661 0.096787 .
freepooryes 0.066277 0.269699 0.246 0.805879
freerepatyes -0.434941 0.166196 -2.617 0.008870 **
nchronicyes 0.109660 0.125380 0.875 0.381779
lchronicyes 0.135612 0.142766 0.950 0.342166
Zero hurdle model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.224621 0.156469 -20.609 < 2e-16 ***
genderfemale 0.305324 0.089648 3.406 0.000660 ***
age 0.700345 0.276089 2.537 0.011191 *
illness1 0.885148 0.136951 6.463 1.02e-10 ***
illness2 1.238227 0.146059 8.478 < 2e-16 ***
illness3 1.263698 0.169344 7.462 8.50e-14 ***
illness4 1.405167 0.195388 7.192 6.40e-13 ***
illness5 1.445585 0.208425 6.936 4.04e-12 ***
reduced 0.154858 0.013488 11.481 < 2e-16 ***
health 0.070464 0.019142 3.681 0.000232 ***
privateyes 0.271192 0.112751 2.405 0.016163 *
freepooryes -0.546177 0.277942 -1.965 0.049406 *
freerepatyes 0.423220 0.153994 2.748 0.005991 **
nchronicyes -0.006256 0.102033 -0.061 0.951106
lchronicyes 0.070658 0.140587 0.503 0.615251
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of iterations in BFGS optimization: 22
Log-likelihood: -2581 on 30 Df
As we see this output has two tables. The above one is for the poisson model performed only on the truncated positive values, and the below one is the result of the logistic regression with only two classes (zero or positive value)
As we did before we plot the results.
predhp<- predict(modelhp,newdata=test[test$visits!=0,],type = "response") plot(test$visits[test$visits!=0],type = "b",col="red") lines(round(predhp),col="blue")
As before by only looking at the plot we can not decide which model is the best. So it is better to use the statistic metrics.
predhp<- predict(modelhp,newdata=test, type = "response") rmsemodelhp<-ModelMetrics::rmse(test$visits,round(predhp)) maemodelhp<-mae(test$visits,round(predhp)) rmsemodelhp [1] 0.7375374 maemodelhp [1] 0.2850242
hurdle model with negative binomial distribution.
Now let’s try to use the negative binomial instead of poisson distribution.
set.seed(123)
modelhnb<-hurdle(visits~.-income, data=train,dist = "negbin")
summary(modelhnb)
Call:
hurdle(formula = visits ~ . - income, data = train, dist = "negbin")
Pearson residuals:
Min 1Q Median 3Q Max
-0.9078 -0.4515 -0.3201 -0.2022 10.6552
Count model coefficients (truncated negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.68462 2.65037 -1.390 0.1645
genderfemale 0.07432 0.17299 0.430 0.6675
age 0.26774 0.54001 0.496 0.6200
illness1 0.27678 0.34564 0.801 0.4233
illness2 0.37093 0.35241 1.053 0.2925
illness3 0.04728 0.39747 0.119 0.9053
illness4 0.40386 0.40517 0.997 0.3189
illness5 0.68213 0.41357 1.649 0.0991 .
reduced 0.15813 0.01935 8.171 3.05e-16 ***
health 0.01891 0.03291 0.575 0.5656
privateyes -0.45711 0.23118 -1.977 0.0480 *
freepooryes 0.03334 0.55282 0.060 0.9519
freerepatyes -0.59189 0.30437 -1.945 0.0518 .
nchronicyes 0.08737 0.21061 0.415 0.6783
lchronicyes 0.30274 0.25846 1.171 0.2415
Log(theta) -2.80552 2.80120 -1.002 0.3166
Zero hurdle model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.224621 0.156469 -20.609 < 2e-16 ***
genderfemale 0.305324 0.089648 3.406 0.000660 ***
age 0.700345 0.276089 2.537 0.011191 *
illness1 0.885148 0.136951 6.463 1.02e-10 ***
illness2 1.238227 0.146059 8.478 < 2e-16 ***
illness3 1.263698 0.169344 7.462 8.50e-14 ***
illness4 1.405167 0.195388 7.192 6.40e-13 ***
illness5 1.445585 0.208425 6.936 4.04e-12 ***
reduced 0.154858 0.013488 11.481 < 2e-16 ***
health 0.070464 0.019142 3.681 0.000232 ***
privateyes 0.271192 0.112751 2.405 0.016163 *
freepooryes -0.546177 0.277942 -1.965 0.049406 *
freerepatyes 0.423220 0.153994 2.748 0.005991 **
nchronicyes -0.006256 0.102033 -0.061 0.951106
lchronicyes 0.070658 0.140587 0.503 0.615251
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta: count = 0.0605
Number of iterations in BFGS optimization: 31
Log-likelihood: -2524 on 31 Df
And let’s plot the difference between the predicted and the actual values of the testing set. .
predhnb<- predict(modelhnb,newdata=test[test$visits!=0,],type = "response") plot(test$visits[test$visits!=0],type = "b",col="red") lines(round(predhnb),col="blue")
And for the metrics.
predhnb<- predict(modelhnb,newdata=test,type = "response") rmsemodelhnb<-ModelMetrics::rmse(test$visits,round(predhnb)) maemodelhnb<-mae(test$visits,round(predhnb)) rmsemodelhnb [1] 0.7408052 maemodelhnb [1] 0.2879227
Zero inflated model
Such as the previous model type , this model also combines two components but with the difference that this model performs a mixture of binomial distribution (between zero and positive values) and the poisson (or negative binomial) distribution for the rest of the values (with the zero included).
Zero inflated model with poisson distribution
Here also we fit tow models one with poisson and one with negative binomial
set.seed(123)
modelzp<-zeroinfl(visits~.-income, data=train,dist = "poisson")
summary(modelzp)
Call:
zeroinfl(formula = visits ~ . - income, data = train, dist = "poisson")
Pearson residuals:
Min 1Q Median 3Q Max
-1.6247 -0.4791 -0.3326 -0.1783 12.3448
Count model coefficients (poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.86466 0.23608 -3.663 0.00025 ***
genderfemale 0.03280 0.09078 0.361 0.71788
age 0.13330 0.26922 0.495 0.62050
illness1 0.32985 0.21846 1.510 0.13107
illness2 0.34799 0.22426 1.552 0.12072
illness3 0.20399 0.24152 0.845 0.39833
illness4 0.44019 0.24324 1.810 0.07034 .
illness5 0.72462 0.23632 3.066 0.00217 **
reduced 0.09679 0.00809 11.964 < 2e-16 ***
health 0.02269 0.01609 1.410 0.15855
privateyes -0.26390 0.12796 -2.062 0.03918 *
freepooryes 0.04860 0.27675 0.176 0.86060
freerepatyes -0.51895 0.17070 -3.040 0.00236 **
nchronicyes 0.08577 0.11489 0.747 0.45536
lchronicyes 0.10876 0.12745 0.853 0.39347
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.44814 0.34932 7.008 2.41e-12 ***
genderfemale -0.48766 0.19038 -2.561 0.010422 *
age -0.88817 0.59030 -1.505 0.132426
illness1 -0.80835 0.31247 -2.587 0.009683 **
illness2 -1.41462 0.35338 -4.003 6.25e-05 ***
illness3 -1.69205 0.44028 -3.843 0.000121 ***
illness4 -1.52225 0.46334 -3.285 0.001018 **
illness5 -1.08745 0.46493 -2.339 0.019338 *
reduced -0.14462 0.03861 -3.746 0.000180 ***
health -0.05795 0.04486 -1.292 0.196418
privateyes -0.73945 0.22597 -3.272 0.001066 **
freepooryes 0.73370 0.41402 1.772 0.076372 .
freerepatyes -1.75457 0.53937 -3.253 0.001142 **
nchronicyes 0.13230 0.22623 0.585 0.558678
lchronicyes 0.03646 0.30619 0.119 0.905217
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of iterations in BFGS optimization: 80
Log-likelihood: -2579 on 30 Df
predzp<- predict(modelzp,newdata=test[test$visits!=0,],type = "response")
plot(test$visits[test$visits!=0],type = "b",col="red")
lines(round(predzp),col="blue")
predzp<- predict(modelzp,newdata=test,type = "response") rmsemodelzp<-ModelMetrics::rmse(test$visits,round(predzp)) maemodelzp<-mae(test$visits,round(predzp)) rmsemodelzp [1] 0.7485897 maemodelzp [1] 0.2898551
Zero inflated model with negative binomial distribution
Let’s this time try the negative binomial distribution.
set.seed(123)
modelznb<-zeroinfl(visits~., data=train,dist = "negbin")
summary(modelznb)
Call:
zeroinfl(formula = visits ~ ., data = train, dist = "negbin")
Pearson residuals:
Min 1Q Median 3Q Max
-1.0440 -0.4582 -0.3031 -0.1680 14.2061
Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.383088 0.241780 -5.720 1.06e-08 ***
genderfemale 0.038982 0.090534 0.431 0.666775
age 0.028836 0.290843 0.099 0.921022
income -0.196594 0.147560 -1.332 0.182763
illness1 0.613261 0.174653 3.511 0.000446 ***
illness2 0.692295 0.179665 3.853 0.000117 ***
illness3 0.664614 0.196063 3.390 0.000699 ***
illness4 0.760167 0.204138 3.724 0.000196 ***
illness5 0.944761 0.206098 4.584 4.56e-06 ***
reduced 0.102651 0.008776 11.697 < 2e-16 ***
health 0.044012 0.015536 2.833 0.004612 **
privateyes -0.168877 0.138680 -1.218 0.223321
freepooryes -0.422737 0.306654 -1.379 0.168034
freerepatyes -0.383568 0.163995 -2.339 0.019341 *
nchronicyes 0.033375 0.107881 0.309 0.757040
lchronicyes 0.065839 0.128987 0.510 0.609748
Log(theta) 0.473940 0.142626 3.323 0.000891 ***
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.280216 0.532170 4.285 1.83e-05 ***
genderfemale -0.726916 0.275312 -2.640 0.00828 **
age -2.002886 0.920168 -2.177 0.02951 *
income -0.180281 0.393267 -0.458 0.64665
illness1 -0.332734 0.347965 -0.956 0.33896
illness2 -1.111993 0.449603 -2.473 0.01339 *
illness3 -0.953342 0.512702 -1.859 0.06296 .
illness4 -1.551433 0.739821 -2.097 0.03599 *
illness5 -1.229820 0.859693 -1.431 0.15256
reduced -1.298180 0.457644 -2.837 0.00456 **
health -0.001444 0.055085 -0.026 0.97909
privateyes -0.817933 0.317762 -2.574 0.01005 *
freepooryes 0.239439 0.664835 0.360 0.71874
freerepatyes -13.563347 519.415149 -0.026 0.97917
nchronicyes 0.045017 0.298234 0.151 0.88002
lchronicyes -0.163689 0.495085 -0.331 0.74093
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta = 1.6063
Number of iterations in BFGS optimization: 68
Log-likelihood: -2512 on 33 Df
predznb<- predict(modelznb,newdata=test,type = "response")
rmsemodelznb<-ModelMetrics::rmse(test$visits,round(predznb))
maemodelznb<-mae(test$visits,round(predznb))
rmsemodelznb
[1] 0.7309579
maemodelznb
[1] 0.2753623
Finally let’s compare all the above models.
rmse<-c(rmsemodelp,rmsemodelqp,rmsemodelnb,rmsemodelhp,rmsemodelhnb,
rmsemodelzp,rmsemodelznb)
mae<-c(maemodelp,maemodelqp,maemodelnb,maemodelhp,maemodelhnb,
maemodelzp,maemodelznb)
models<-c("pois","q_pois","nb","h_pois","h_nb","zer_pois","zer_nb")
data.frame(models,rmse,mae)%>%
arrange(rmse)
models rmse mae
1 zer_nb 0.7309579 0.2753623
2 h_pois 0.7375374 0.2850242
3 pois 0.7381921 0.2840580
4 q_pois 0.7381921 0.2840580
5 h_nb 0.7408052 0.2879227
6 zer_pois 0.7485897 0.2898551
7 nb 0.7808085 0.2966184
Both metrics have chosen the zero inflated negative binomial model as the best model with minimum rmse value 0.7309579, and minimum mae value 0.2753623. this result is in line with the fact that this kind of models take care of the zero inflated data and at the same time the overdispersion problem.
Conclusion:
If the data is truly follows Poisson distribution then all the the other models have extra parameters that, during training process, converges to the optimum parameter values for poisson, this relation is like the linear regression to the Generealized least squares. However, if the data is very skewed towards zero then it should be better to use the last two models to take care of this issue.
Furhter reading:
- Michael J. Crawley, The R book, WILEY, UK, 2013. http://www.bio.ic.ac.uk/research/mjcraw/therbook/index.htm
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.