Change in consumer surplus in R

I want to share a short code to calculate change in consumer surplus in R after a regression. I will borrow the last blog post’s example to show how to do it.

My goal is to understand how consumer surplus (CS) changes when a fishing port is present. Calculating change in consumer surplus is not as straightforward as calculating the change in CS, but still quite simple.

The equation to calculate the change in consumer surplus is written in Hanley et al. (2003). They wanted to calculate change in consumer surplus due to coastal water quality improvements. Basically, the change in consumer surplus is the change of the expected number of trips with and without the quality change (such as presence of a fishing port). In my application on fishing presence problem, the change in consumer surplus is:

 \Delta CS = \frac{exp(intercept-beta_{TC}*TC + beta_X*X +beta_{Fishpres} * Fishpres )}{-beta_{TC}} - \frac{exp(intercept-beta_{TC}*TC + beta_X*X +beta_{Fishpres} * 0)}{-beta_{TC}}

where X denotes a matrix of relevant explanatory variables. The first numerator denotes the expected number of trips with fishing presence and the second numerator is the expected number of trips without fishing presence.

Let me assume the matrix of explanatory variables include only substitute travel cost and income. I need the mean value to fill in the change in consumer surplus equation:

> mean(Dataset$TC, na.rm = TRUE)
[1] 6.882093
> mean(Dataset$substitute_TC, na.rm = TRUE)
[1] 73.51123
> mean(Dataset$income, na.rm = TRUE)
[1] 19236.86

I can substitute the average values of travel cost (TC), Substitute Travel cost (SubsTC) and income in the equation to find the change in consumer surplus from fishing presence for the average individual:

 \Delta CS = \frac{exp(intercept-beta_{TC}*6.882 + beta_X*73.5112+beta_{income}*19237 +beta_{Fishpres} * 1 )}{-beta_{TC}} - \frac{exp(intercept-beta_{TC}*6.882 + beta_X*73.5112+beta_{income}*19237)}{-beta_{TC}}

All we are missing are the parameter values. To obtain them, I run a truncated negative binomial model on my count data using the VGAM package:

> m1 <- vglm(visits ~ TC + 
+                substitute_TC +
+                income +
+                fishing_presence 
+            , family = posnegbinomial(), data = Dataset)
There were 50 or more warnings (use warnings() to see the first 50)
> summary(m1)

                      Estimate Std. Error z value Pr(>|z|)    
(Intercept):1       -7.515e-01  5.140e-01  -1.462   0.1437    
(Intercept):2       -2.833e+00  5.651e-01  -5.014 5.34e-07 ***
TC                  -3.313e-02  3.638e-03  -9.108  < 2e-16 ***
substitute_TC        4.568e-04  2.118e-04   2.157   0.0310 *  
income               9.289e-06  3.278e-06   2.834   0.0046 ** 
fishing_presence     1.681e-01  9.993e-02   1.682   0.0926 .  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The coefficients from a vglm object are stored in the m1@coefficients object. Using these coefficients, if I calculate the consumer surplus by hand, I get:

 \Delta CS = \frac{exp(-7.515302e-01 - 6.882093*3.312949e-02 + 73.51123*4.568340e-04 + 19236.86*9.289449e-06 + 1*1.680805e-01)}{3.312949e-02} - \frac{exp(-7.515302e-01 - 6.882093*3.312949e-02 + 73.51123*4.568340e-04 + 19236.86*9.289449e-06)}{3.312949e-02} = 2.565

This would imply that the change in consumer surplus is around 2.56 euros from presence of a fishing port per visit.

It is worth exploring how to automate this equation in case I add more explanatory variables, or if I want to look into a different variable of interest. Therefore, I am sharing below a code to calculate consumer surplus after running a truncated negative binomial model:

### Change in consumer surplus function

change_CS<-function(df,model,tc,X, level_0, level_1){
  means <- 0
  for(i in 3:length(coef(model))){
       means[i] <- mean(df[[names(coef(model))[i]]], na.rm=TRUE)
  names(means) <- names(coef(m1))
  means[] <- 0
  means[X] <- 0
  delta_CS <-  exp( coef(model)[1] + sum( means*coef(model) ) + level_1*coef(model)[grep(X, names(coef(model)))])/
               (-coef(model)[grep(tc, names(coef(model)))] )   -
               exp( coef(model)[1] + sum( means*coef(model) ) + level_0*coef(model)[grep(X, names(coef(model)))])/
               (-coef(model)[grep(tc, names(coef(model)))] )           

The researcher needs to define six arguments: 1) the dataset with the estimated data; 2) the model name, 3) the cost variable name, 4) the name of the variable of interest, and 5) the value of the variable of interest before and 6) after the change.

In my case, I run the following command to run this function:

> change_CS(Dataset,m1,"\\bTC\\b","fishing_presence", 0, 1)

As expected, I obtain the exact same change in CS value as doing it by hand. I have confirmed that this function works if I expand the matrix of explanatory variables, as well as if I change the name of the variable of interest. This function only works for the truncated NB model (for now), but it can easily be adapted for other count data models.


Hanley, N., Bell, D., & Alvarez-Farizo, B. (2003). Valuing the benefits of coastal water quality improvements using contingent and real behaviour. Environmental and resource economics24(3), 273-285.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s