Economic Experiments for Non-Economists Using R
Clinton L. Neill, Ph.D.
7/10/2020
Code for Economic Experiments
Within this code you will find a step by step guide of completing an economic experiment to measure willingness-to-pay (WTP). This document is to serve as a guide to the basics and not an all inclusive document. Note that the packages and code presented here are generally flexible enough to perform more advanced analysis. However, some advanced models may require you to utilize other programs such as MATLAB, SAS, or LIMDEP. These other programs are also better at handling extremely large data sets. I am happy to provide code and guidance if you need to pursue other programs.
Basics about how the document is formatted and presented
To give you some basics about the document, the following is an example of R code and output.
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Note that the output is denoted by the ## on each line.
Also the text style is also different from regular text. I have used R
Markdown to create this document to easily introduce code and output to
you throughout this tutorial.
Including Plots
Sometimes I will also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.
Note that I sometimes suppress output for easier viewing of the document. If you have the correct data and code, you should be able to run the code and see the output yourself.
Defining your Question
The first step in any economic experiment it to define your question. What economic phenomenon are you interested in studying? Why is it important that you study this question? Who/What will it impact? I will work through an example of an economic question throughout this document, but I encourage you to think of your own economic question and code with that in mind.
Remember that economics is the study of trade offs, choices, and scarcity. Most of the choices we make everyday have an economic impact on ourselves, our local economies, and beyond. By defining your question before you even think about the experimental design will help facilitate more robust research.
The Science of the Problem
Once you have defined your economic problem, the first step is to have a scientific way to design your experiment. This helps to ensure that we are creating unbiased tests for our hypothesis. Without this vital step of the process we would bias our own results and our economic recommendations would be misleading. This is perhaps the most pronounced flaw I, personally, find when reading WTP studies performed by non-economists. Simply asking whether someone is willing-to-pay for a product invites response bias and usually inflates WTP values to non-nonsensical levels. Economists have developed many methods to reduce these issues and helps us to have scientifically backed conclusions.
Sample Size Calculations
The first step to experimental design is to determine your sample size. Depending on your economic question, the population that is affected by your economic question will vary. For example, if we are interested in how a change in U.S. food safety regulations affects the WTP for fresh produce, the population of interest is everyone in the U.S. (or >325 million people). On the other hand, if we are interested in how WTP for a new craft beer, that population is much smaller (~40% of those who are of the legal drinking age - https://www.statista.com/statistics/957412/craft-beer-drinker-share-us/). Populations can be even smaller if you are concerned with people in a single state, county, or city/town.
Statistics has provided us with equations to calculate how many people would need to be surveyed to be confident in that the behavior we observe from survey responses is representative of the entire population. Note that these calculations have asymptotic properties - meaning that once the population of interest reaches a certain level, the sample size needed remains the same. The common knowledge is that 1000 is sufficient to measure a phenomenon within a 95%-99% confidence level for sample sizes of 100,000 people or more.
While I have not found a package that calculates sample size for you, I have written a function in R to calculate it for you:
sample.size.table = function(margin=.5, c.interval=.05, population) {
z.val=c(1.281551565545, 1.644853626951, 1.959963984540,
2.326347874041, 2.575829303549, 2.807033768344,
3.090232306168, 3.290526731492, 3.890591886413)
ss = (z.val^2 * margin * (1-margin))/(c.interval^2)
p.ss = ss/(1 + ((ss-1)/population))
c.level = c("80%","90%","95%","98%","99%",
"99.5%","99.8%","99.9%","99.99%")
results = data.frame(c.level, round(p.ss, digits = 0))
names(results) = c("Confidence Level", "Sample Size")
METHOD = c("Suggested sample sizes at different confidence levels")
moe = paste((c.interval*100), "%", sep="")
resp.dist = paste((margin*100),"%", sep="")
pre = structure(list(Population=population,
"Margin of error" = moe,
"Response distribution" = resp.dist,
method = METHOD),
class = "power.htest")
print(pre)
print(results)
}All you need to do is input your response margin (the percentage of the population you think will respond to your treatment - if unknown leave at 50%), the margin or error (how much random sampling error - the default is 5%), and the population. If I want to leave my response margin and margin of error at the default, then the sample size for a population of 1,000 is:
sample.size.table(, , 1000)##
## Suggested sample sizes at different confidence levels
##
## Population = 1000
## Margin of error = 5%
## Response distribution = 50%
##
## Confidence Level Sample Size
## 1 80% 141
## 2 90% 213
## 3 95% 278
## 4 98% 351
## 5 99% 399
## 6 99.5% 441
## 7 99.8% 489
## 8 99.9% 520
## 9 99.99% 602As you can see, the R output returns a table to show how your chosen confidence interval (95%) compares to other confidence levels. You would need to survey 278 people to observe your phenomenon with 95% confidence if you population was 1,000 people. Sampling more than this number will improve your analysis results, but with diminishing marginal returns to accuracy - unless your assumed response distribution is incorrect. I would air on the side of caution if you don’t know your response distribution and over sample if able.
Experimental Designs Common for Economics
The most common types of experiments for economics are choice-based conjoint (also referred to more generally as choice experiments), contingent valuation, and best-worst experiments. Choice-based conjoint and best-worst experiments require you to use a statistical experimental design to ensure that the treatments are not correlated and can be identified (rather than confounded or aliased).
Choice Experiments
The most common way to design choice experiments is to a full or
fractional factorial. This type of design allows for design
orthogonality (no correlation) between treatments and offers a degree of
randomization between the choice options. There are a couple packages
you can use to render the designs: AlgDesign and
conjoint. The former allows for more flexibility and
calculates design efficiency for you. However, designs that require
higher numbers of levels for the treatments (more than 2 or 3), a lot of
manual manipulation is required. The conjoint package is
much easier to use and can handle treatments with greater numbers of
levels.
Some quick definitions/notations: treatments are the general things
you want to test and levels are the variations within treatments. For
example, if I wanted to test WTP for milk fat content I would have two
treatments: milk fat and price. Within the milk fat treatment I have 3
levels - skim, 2%, and whole; similarly, I can also specify 3 levels of
the price treatment - 2.99, 3.99, and 4.99. Let’s use this example to
create a choice experiment using the conjoint package.
#install.packages("conjoint")
library(conjoint)## Warning: package 'conjoint' was built under R version 4.2.3Assuming we only have one choice option, the one we design or a ‘no buy’ option, we can determine the number of options a full factorial would contain. Assume L=“levels” and T=“Treatments”, then the formula for the full factorial is L^(# of T with same # of L). Since we have two treatments each with 3 levels: 3^2 = 9 total choice questions. However, if we reduced the number of price levels to 2, then the math is a little different: 3^(1)*2(1) = 6 total choices. Let’s use our original 2 treatments both at 3 levels to create a full factorial:
We use expand.grid to allow us to create an easy
data.frame for our attributes and levels
experiment<-expand.grid( price=c("$2.99","$3.99","$4.99"),
FatContent=c("Skim","2%","Whole"))Now to create the design we use the
caFactorialDesign(data, type='', seed=' ') to create the
design. The type is to specify a full or fractional
factorial, while seed tells you the random set to start the
algorithm. Let’s start with a full factorial.
design=caFactorialDesign(data=experiment,type="full", seed=55751)
print(design) ## price FatContent
## 1 $2.99 Skim
## 2 $3.99 Skim
## 3 $4.99 Skim
## 4 $2.99 2%
## 5 $3.99 2%
## 6 $4.99 2%
## 7 $2.99 Whole
## 8 $3.99 Whole
## 9 $4.99 Wholeprint(cor(caEncodedDesign(design))) ## price FatContent
## price 1 0
## FatContent 0 1As we can see, all 9 combinations are presented. The
cor(caEncodedDesign(design)) just tells us whether their
correlation between treatments. If the correlation is high, you likely
need a different design as you might be confounding/aliasing treatment
effects.
However, this would only be one ‘choice option’ and most choice experiments have at least two choice options in addition to a ‘no buy’ option. To expand this to multiple ‘choice options’, simply double the attributes in the experiment:
experiment2<-expand.grid( price1=c("$2.99","$3.99","$4.99"),
FatContent1=c("Skim","2%","Whole"),
price2=c("$2.99","$3.99","$4.99"),
FatContent2=c("Skim","2%","Whole"))
design2=caFactorialDesign(data=experiment2,type="full")
print(design2) ## price1 FatContent1 price2 FatContent2
## 1 $2.99 Skim $2.99 Skim
## 2 $3.99 Skim $2.99 Skim
## 3 $4.99 Skim $2.99 Skim
## 4 $2.99 2% $2.99 Skim
## 5 $3.99 2% $2.99 Skim
## 6 $4.99 2% $2.99 Skim
## 7 $2.99 Whole $2.99 Skim
## 8 $3.99 Whole $2.99 Skim
## 9 $4.99 Whole $2.99 Skim
## 10 $2.99 Skim $3.99 Skim
## 11 $3.99 Skim $3.99 Skim
## 12 $4.99 Skim $3.99 Skim
## 13 $2.99 2% $3.99 Skim
## 14 $3.99 2% $3.99 Skim
## 15 $4.99 2% $3.99 Skim
## 16 $2.99 Whole $3.99 Skim
## 17 $3.99 Whole $3.99 Skim
## 18 $4.99 Whole $3.99 Skim
## 19 $2.99 Skim $4.99 Skim
## 20 $3.99 Skim $4.99 Skim
## 21 $4.99 Skim $4.99 Skim
## 22 $2.99 2% $4.99 Skim
## 23 $3.99 2% $4.99 Skim
## 24 $4.99 2% $4.99 Skim
## 25 $2.99 Whole $4.99 Skim
## 26 $3.99 Whole $4.99 Skim
## 27 $4.99 Whole $4.99 Skim
## 28 $2.99 Skim $2.99 2%
## 29 $3.99 Skim $2.99 2%
## 30 $4.99 Skim $2.99 2%
## 31 $2.99 2% $2.99 2%
## 32 $3.99 2% $2.99 2%
## 33 $4.99 2% $2.99 2%
## 34 $2.99 Whole $2.99 2%
## 35 $3.99 Whole $2.99 2%
## 36 $4.99 Whole $2.99 2%
## 37 $2.99 Skim $3.99 2%
## 38 $3.99 Skim $3.99 2%
## 39 $4.99 Skim $3.99 2%
## 40 $2.99 2% $3.99 2%
## 41 $3.99 2% $3.99 2%
## 42 $4.99 2% $3.99 2%
## 43 $2.99 Whole $3.99 2%
## 44 $3.99 Whole $3.99 2%
## 45 $4.99 Whole $3.99 2%
## 46 $2.99 Skim $4.99 2%
## 47 $3.99 Skim $4.99 2%
## 48 $4.99 Skim $4.99 2%
## 49 $2.99 2% $4.99 2%
## 50 $3.99 2% $4.99 2%
## 51 $4.99 2% $4.99 2%
## 52 $2.99 Whole $4.99 2%
## 53 $3.99 Whole $4.99 2%
## 54 $4.99 Whole $4.99 2%
## 55 $2.99 Skim $2.99 Whole
## 56 $3.99 Skim $2.99 Whole
## 57 $4.99 Skim $2.99 Whole
## 58 $2.99 2% $2.99 Whole
## 59 $3.99 2% $2.99 Whole
## 60 $4.99 2% $2.99 Whole
## 61 $2.99 Whole $2.99 Whole
## 62 $3.99 Whole $2.99 Whole
## 63 $4.99 Whole $2.99 Whole
## 64 $2.99 Skim $3.99 Whole
## 65 $3.99 Skim $3.99 Whole
## 66 $4.99 Skim $3.99 Whole
## 67 $2.99 2% $3.99 Whole
## 68 $3.99 2% $3.99 Whole
## 69 $4.99 2% $3.99 Whole
## 70 $2.99 Whole $3.99 Whole
## 71 $3.99 Whole $3.99 Whole
## 72 $4.99 Whole $3.99 Whole
## 73 $2.99 Skim $4.99 Whole
## 74 $3.99 Skim $4.99 Whole
## 75 $4.99 Skim $4.99 Whole
## 76 $2.99 2% $4.99 Whole
## 77 $3.99 2% $4.99 Whole
## 78 $4.99 2% $4.99 Whole
## 79 $2.99 Whole $4.99 Whole
## 80 $3.99 Whole $4.99 Whole
## 81 $4.99 Whole $4.99 Wholeprint(cor(caEncodedDesign(design2)))## price1 FatContent1 price2 FatContent2
## price1 1 0 0 0
## FatContent1 0 1 0 0
## price2 0 0 1 0
## FatContent2 0 0 0 1As you can see, this rapidly expands the total number of total options in the design. This is obviously quite cumbersome for any one person to evaluate. In addition, there are scenarios that have the same choice option for both options - see the 81st choice scenario:
design2[81,]This doesn’t provide us any information about one’s choices and we often eliminate choice questions with the same choice twice. To reduce experimental design size and remain efficient let’s use a fractional factorial design:
design3.fract=caFactorialDesign(data=experiment2,type="fractional")
print(design3.fract)## price1 FatContent1 price2 FatContent2
## 2 $3.99 Skim $2.99 Skim
## 15 $4.99 2% $3.99 Skim
## 16 $2.99 Whole $3.99 Skim
## 24 $4.99 2% $4.99 Skim
## 26 $3.99 Whole $4.99 Skim
## 30 $4.99 Skim $2.99 2%
## 35 $3.99 Whole $2.99 2%
## 41 $3.99 2% $3.99 2%
## 45 $4.99 Whole $3.99 2%
## 46 $2.99 Skim $4.99 2%
## 50 $3.99 2% $4.99 2%
## 58 $2.99 2% $2.99 Whole
## 65 $3.99 Skim $3.99 Whole
## 81 $4.99 Whole $4.99 Wholeprint(cor(caEncodedDesign(design3.fract)))## price1 FatContent1 price2 FatContent2
## price1 1.00000000 0.1032796 0.1032796 -0.09259259
## FatContent1 0.10327956 1.0000000 0.2160000 -0.10327956
## price2 0.10327956 0.2160000 1.0000000 -0.10327956
## FatContent2 -0.09259259 -0.1032796 -0.1032796 1.00000000As you can see the design only returns specific rows of the full factorial. Again, duplicate choice options in the same scenario (row) are not useful and can be eliminated at this stage or earlier (earlier elimination require more coding that is not presented here). One thing to note, you need to ensure that all of the treatment levels are present across the fractional design. If all treatment levels are not present, you are missing variation and we prefer at least one direct observation of that level within the choice experiment.
As you can see the cross-correlations are very small which means this design does not have higher order effects.
In this example, 2 rows have duplicate choices (41 and 81). If we
change the type='' value the optimal design will
change:
design3.1.fract=caFactorialDesign(data=experiment2,type="null")
print(design3.1.fract)## price1 FatContent1 price2 FatContent2
## 6 $4.99 2% $2.99 Skim
## 16 $2.99 Whole $3.99 Skim
## 20 $3.99 Skim $4.99 Skim
## 29 $3.99 Skim $2.99 2%
## 36 $4.99 Whole $2.99 2%
## 39 $4.99 Skim $3.99 2%
## 40 $2.99 2% $3.99 2%
## 50 $3.99 2% $4.99 2%
## 52 $2.99 Whole $4.99 2%
## 55 $2.99 Skim $2.99 Whole
## 59 $3.99 2% $2.99 Whole
## 71 $3.99 Whole $3.99 Whole
## 75 $4.99 Skim $4.99 Whole
## 78 $4.99 2% $4.99 Wholeprint(cor(caEncodedDesign(design3.1.fract)))## price1 FatContent1 price2 FatContent2
## price1 1.0000000 -0.2160000 0 0.1032796
## FatContent1 -0.2160000 1.0000000 0 -0.1032796
## price2 0.0000000 0.0000000 1 0.0000000
## FatContent2 0.1032796 -0.1032796 0 1.0000000By changing type to ‘orthogonal’, we obtain a design that has zero cross attribute correlation. This is typically what you want, but not necessary if cross correlations are small.
design3.ortho.fract=caFactorialDesign(data=experiment2,type="orthogonal", seed=68795)
print(design3.ortho.fract)## price1 FatContent1 price2 FatContent2
## 5 $3.99 2% $2.99 Skim
## 18 $4.99 Whole $3.99 Skim
## 19 $2.99 Skim $4.99 Skim
## 30 $4.99 Skim $2.99 2%
## 40 $2.99 2% $3.99 2%
## 53 $3.99 Whole $4.99 2%
## 61 $2.99 Whole $2.99 Whole
## 65 $3.99 Skim $3.99 Whole
## 78 $4.99 2% $4.99 Wholeprint(cor(caEncodedDesign(design3.ortho.fract)))## price1 FatContent1 price2 FatContent2
## price1 1 0 0 0
## FatContent1 0 1 0 0
## price2 0 0 1 0
## FatContent2 0 0 0 1Once you have the statistical design created and duplicates removed, you can design your experiment in Qualtrics or other program. This is where the “Art” of the process comes in. There are plenty of samples of visually appealing surveys on the internet. If you would like to see one of mine, just ask!
Best-Worst Experiments
An interesting type of economic experiment is a Best-Worst study that
forces survey participants to make trade offs between
alternatives/treatments. This type of experiment helps in creating
rankings and does not allow people to avoid making a choice (choice
experiments typically allow for a ‘no buy’ option). To create a
Best-Worst study we use a Balanced Incomplete Block Design (BIBD). For
this type of design we take advantage of the crossdes
package:
#install.packages("crossdes")
library(crossdes)## Warning: package 'crossdes' was built under R version 4.2.3## Loading required package: AlgDesign## Loading required package: gtools## Warning: package 'gtools' was built under R version 4.2.3Without solving some combonatorics math, BIBDs are difficult to
create appropriately. See the following document (specifically page 13-1
- which is the second page) https://www.stat.purdue.edu/~bacraig/notes1/topic13.pdf
for the specific equations. The function we use for BIBDs is the
find.BIB('treatments', 'blocks', 'elements). Without using
the combonatorics math tricks, we can attempt to create a BIBD through
trial and error as follows:
find.BIB(5, 4, 3) ## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 1 3 4
## [3,] 2 4 5
## [4,] 1 2 3isGYD(find.BIB(5, 4, 3), tables = TRUE, type=TRUE) ##
## [1] The design is neither balanced w.r.t. rows nor w.r.t. columns.
##
## $`Number of occurrences of treatments in d`
## 1 2 3 4 5
## 3 2 2 2 3
##
## $`Row incidence matrix of d`
## 1 2 3 4
## 1 1 1 0 1
## 2 1 0 1 0
## 3 0 0 1 1
## 4 0 1 1 0
## 5 1 1 0 1
##
## $`Column incidence matrix of d`
## 1 2 3
## 1 3 0 0
## 2 1 1 0
## 3 0 2 0
## 4 0 1 1
## 5 0 0 3
##
## $`Concurrence w.r.t. rows`
## 1 2 3 4 5
## 1 3 1 1 1 3
## 2 1 2 1 1 1
## 3 1 1 2 1 1
## 4 1 1 1 2 1
## 5 3 1 1 1 3
##
## $`Concurrence w.r.t. columns`
## 1 2 3 4 5
## 1 9 3 0 0 0
## 2 3 2 2 1 0
## 3 0 2 4 2 0
## 4 0 1 2 2 3
## 5 0 0 0 3 9We attempted to create the following design: five treatments in four
blocks of three elements. The isGYD() function tests to see
if the design is Balanced. Obviously this is not a balanced design. If
you think of each row as a block there is clearly not a balance between
the treatments. If we had used the combonatorics math, we would have
found that we need 10 blocks to create a balanced design:
find.BIB(5, 10, 3) ## [,1] [,2] [,3]
## [1,] 1 4 5
## [2,] 2 3 5
## [3,] 1 2 4
## [4,] 3 4 5
## [5,] 2 3 4
## [6,] 1 2 5
## [7,] 1 3 5
## [8,] 1 2 3
## [9,] 2 4 5
## [10,] 1 3 4isGYD(find.BIB(5, 10, 3), tables = TRUE, type=TRUE) ##
## [1] The design is a balanced incomplete block design w.r.t. rows.
##
## $`Number of occurrences of treatments in d`
## 1 2 3 4 5
## 6 6 6 6 6
##
## $`Row incidence matrix of d`
## 1 2 3 4 5 6 7 8 9 10
## 1 1 1 1 0 1 0 1 0 0 1
## 2 0 0 1 0 1 1 1 1 1 0
## 3 1 0 1 1 0 1 0 0 1 1
## 4 1 1 0 1 0 0 1 1 1 0
## 5 0 1 0 1 1 1 0 1 0 1
##
## $`Column incidence matrix of d`
## 1 2 3
## 1 6 0 0
## 2 3 3 0
## 3 1 4 1
## 4 0 3 3
## 5 0 0 6
##
## $`Concurrence w.r.t. rows`
## 1 2 3 4 5
## 1 6 3 3 3 3
## 2 3 6 3 3 3
## 3 3 3 6 3 3
## 4 3 3 3 6 3
## 5 3 3 3 3 6
##
## $`Concurrence w.r.t. columns`
## 1 2 3 4 5
## 1 36 18 6 0 0
## 2 18 18 15 9 0
## 3 6 15 18 15 6
## 4 0 9 15 18 18
## 5 0 0 6 18 36This design would be considered a BIBD.
A Side Note about BIBDs and Choice Experiments
If you think about each potential combination of all attributes and levels, you can create a choice experiment with a BIBD instead of a fractional factorial design. However, if the design is large enough it would be impossible to get a fully balanced design. Instead, we could use a Partially Balanced Incomplete Block Design (PBIBD), which would require every option to appear an equal number of times throughout the choice experiment - though the occurrence of every pair does not appear the same number of times as in a BIBD.
No matter the design you end up using, it is important to save the final output. You will need this output when you do analysis.
Exercise
To prepare for the next section, your goal is now to create a choice experiment using a factorial design. This is based on a real experiment that was done last year. We were interested in how consumers valued production and sensory attributes in new varieties of edamame. Production was varied at 2 levels (Non-GMO & Organic); Sensory attributes was varied at 5 levels (Sweet, Salty, Bitter, Nutty, Grassy); and Price at 3 levels (3.25, 4.25, and 5.25 - this was in terms of price for a salad that contained edamame).
- Without using R to create the design, how many choice questions would be in a full factorial if there are 2 choice options in each choice question?
- Generate the full factorial
- Generate the fractional factorial
- How many choice questions does this design suggest using in our experiment?
Answering your Economic Question
Now that you have determined your sample size, designed a proper experiment, and collected your data it is finally time to analyze your data. This is usually the most exciting part of the process because you can finally see all your hard work pay off. You can now answer the question you develop at the beginning.
The statistical analysis of your data is heavily dependent on your experimental design. Often times we think about what statistical model we want to run before we do any experimental design. It is a bit backwards but necessary as you need the proper data to run the proper model. Within this section I will cover, in detail, how to use a multinomial/conditional logit, a random parameters/mixed logit, and a censored regression. I will also discuss the latent class model for choice experiments, but do not provide code. The various logit models are explicitly used for choice experiments and the censored regression is used for contingent valuation experiments. To aid you in learning the models, I will provide data to use in the analysis. The fractional factorial design you created for the edamame exercise is similar to the one used for the edamame choice experiment data you will analyze in this section. A separate edamame WTP experiment was conducted using the contingent valuation method and will that data will be used for the censored regression example. More detail on each experiment is available in those sections below.
Data Manipulation
As with all survey data that is collected via a platform, the data is usually never in the format we need for analysis. Choice experiment data in particular requires specific formatting, and this various by statistical program. Here I will present the data manipulations needed to be done to analyze the choice data in R. Note that if you need to analyze the data in a different software, I can provide code for that as well. Also, the data for the censored regression was collected via paper survey and was formatted during data input in the correct format.
As always, ensure that your working directory is set to where your data is stored for faster data importation. You can check your working directory as follows:
getwd()## [1] "C:/Users/cneill/Dropbox/Crash Course - Econ Exp"If the working directory is not where your data is stored you can change it using the following code (note that this is where I have the data stored and I am using windows):
setwd("C:\\Users\\cneill\\Dropbox\\Crash Course - Econ Exp")The data I provided you is in an excel document so you will need to
install the readxl package. Also, a couple functions used
are from the tidyr package.
#install.packages(c("realxl", "tidyr"))
library(readxl)
library(tidyr)## Warning: package 'tidyr' was built under R version 4.2.3Now we can import the data. We need to import two data files: (1) the attributes file which is the experimental design you created in from the last section; and (2) the actual survey data from qualtrics:
edamame_attributes <- read_excel("Attributes_Edamame_mlogit_r.xlsx")
edamame_survey <- read_excel("Edamame Example.xlsx")Note that the data is in a ‘wide’ format. The data needs to be transformed to be in a ‘long’ format.
edamame_survey edamame_long <- gather(edamame_survey, Q_num, decision,
Q1:Q15, factor_key=TRUE)
edamame_longNow we need to reorganize the data to be by person by question rather than by question by person.
edamame_long <- edamame_long[order(edamame_long$id),]Next, we create a choiceid variable in
edamame_long so that we can merge the two data sets
choiceid <- substring(edamame_long$Q_num, 2)
edamame_long$choiceid <- choiceidThird, we merge the two data frames and reorder the data frame again.
edamame_ce <- merge(edamame_attributes, edamame_long, by = 'choiceid')
edamame_ce <- edamame_ce[order(edamame_ce$id, edamame_ce$choiceid),]Finally, we need to identify the ‘choice’ made in each scenario.
edamame_ce$choice <- ifelse(edamame_ce$decision == 1, "A",
ifelse(edamame_ce$decision ==2, "B", "C"))Huzzah! We now have the data in the format we need to do analysis!
Data Analysis
As I mentioned before, now we can analyze our data. It is important that you review the math and theory behind the models to ensure you know what the model is doing. If you need further explanation, reach out so that I can ensure that the computer is not acting as a ‘black box’. Knowing the math and theory will ensure you are properly interpreting the results.
In the following sections, you will need the following packages:
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("Icens") ## package 'BiocVersion' successfully unpacked and MD5 sums checked
## package 'Icens' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\cneill\AppData\Local\Temp\1\RtmpGswpdA\downloaded_packages### If you are installing this for the first time, you will need to enter the letter "a" into the console. This is being pulled from outside of CRAN and is necessary for the mlogit() function to work.
#install.packages(c("mlogit", "DCchoice", "dfidx", "realxl", "interval"))
library(DCchoice)## Warning: package 'DCchoice' was built under R version 4.2.3library(interval)## Warning: package 'interval' was built under R version 4.2.3## Warning: package 'perm' was built under R version 4.2.3library(mlogit)## Warning: package 'mlogit' was built under R version 4.2.3## Warning: package 'dfidx' was built under R version 4.2.3library(dfidx)Data exploration
Looking at some basic statistics of your data is always the ideal place to start before any actual analysis is done. It gives us the chance to form hypothesis about model parameters so no surprises come up.
Let’s look at some basic statistics first:
summary(edamame_ce)## choiceid price_A sweet_A salty_A bitter_A
## Min. : 1 Min. :3.25 Min. :0.0 Min. :0.0 Min. :0.0
## 1st Qu.: 4 1st Qu.:3.25 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.0
## Median : 8 Median :4.25 Median :0.0 Median :0.0 Median :0.0
## Mean : 8 Mean :4.25 Mean :0.2 Mean :0.2 Mean :0.2
## 3rd Qu.:12 3rd Qu.:5.25 3rd Qu.:0.0 3rd Qu.:0.0 3rd Qu.:0.0
## Max. :15 Max. :5.25 Max. :1.0 Max. :1.0 Max. :1.0
##
## nutty_A grassy_A non-gmo_A organic_A price_B
## Min. :0.0 Min. :0.0 Min. :0.0000 Min. :0.0000 Min. :3.25
## 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:3.25
## Median :0.0 Median :0.0 Median :1.0000 Median :0.0000 Median :4.25
## Mean :0.2 Mean :0.2 Mean :0.5333 Mean :0.4667 Mean :4.25
## 3rd Qu.:0.0 3rd Qu.:0.0 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:5.25
## Max. :1.0 Max. :1.0 Max. :1.0000 Max. :1.0000 Max. :5.25
##
## sweet_B salty_B bitter_B nutty_B grassy_B
## Min. :0.0 Min. :0.0 Min. :0.0 Min. :0.0 Min. :0.0
## 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.0 1st Qu.:0.0
## Median :0.0 Median :0.0 Median :0.0 Median :0.0 Median :0.0
## Mean :0.2 Mean :0.2 Mean :0.2 Mean :0.2 Mean :0.2
## 3rd Qu.:0.0 3rd Qu.:0.0 3rd Qu.:0.0 3rd Qu.:0.0 3rd Qu.:0.0
## Max. :1.0 Max. :1.0 Max. :1.0 Max. :1.0 Max. :1.0
##
## non-gmo_B organic_B price_C sweet_C salty_C
## Min. :0.0000 Min. :0.0000 Min. :0 Min. :0 Min. :0
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :1.0000 Median :0.0000 Median :0 Median :0 Median :0
## Mean :0.6667 Mean :0.3333 Mean :0 Mean :0 Mean :0
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :1.0000 Max. :1.0000 Max. :0 Max. :0 Max. :0
##
## bitter_C nutty_C grassy_C non-gmo_C organic_C id
## Min. :0 Min. :0 Min. :0 Min. :0 Min. :0 Min. : 1
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.: 41
## Median :0 Median :0 Median :0 Median :0 Median :0 Median : 82
## Mean :0 Mean :0 Mean :0 Mean :0 Mean :0 Mean : 82
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:123
## Max. :0 Max. :0 Max. :0 Max. :0 Max. :0 Max. :163
##
## Q_num decision choice
## Q1 : 163 Min. :1.000 Length:2445
## Q2 : 163 1st Qu.:1.000 Class :character
## Q3 : 163 Median :2.000 Mode :character
## Q4 : 163 Mean :1.657
## Q5 : 163 3rd Qu.:2.000
## Q6 : 163 Max. :3.000
## (Other):1467Multinomial/Conditional Logit (MNL/CL)
The multinomial/conditional logit model is always the first model you should run when performing choice experiment analysis. This model is the foundation of the more complex model and gives us a pretty good indication on what type of advanced model might fit our data better.
Before running the choice model we need to perform a couple more data
transformations that are specific to running the model with the
mlogit package. Creating choiceid2 is
necessary as the model analyzes every choice made individually and then
looks across each of those choice by person:
edamame_ce$choiceid2 <- 1:nrow(edamame_ce)
edamame.mlogit <- dfidx(edamame_ce, choice = "choice", varying = 2:25,
sep = "_", idx = list(c("choiceid2", "id")), idnames = c("chid", "alt"))Now we can perform the conditional logit analysis use the
mlogit() function. Remember that all of the variables are
coded as dummy variables, so much drop one of sensory attributes (sweet)
in order for the model to run. We do not need to drop one of the
production attributes (non-gmo or organic) as this is relative to the
‘no buy’ option that was present in the survey as was the third (C)
option in the data (attributes) file. You could reverse which level you
drop from either of the variables/treatments, just remember this will
change your interpretation. You should drop a level from whichever
treatment makes the most sense. Estimation of the model is:
mnl.edamame <- mlogit(choice ~ 0 + `non-gmo` + organic + salty +bitter + nutty + grassy + price,
data=edamame.mlogit)Now to look at the results
summary(mnl.edamame)##
## Call:
## mlogit(formula = choice ~ 0 + `non-gmo` + organic + salty + bitter +
## nutty + grassy + price, data = edamame.mlogit, method = "nr")
##
## Frequencies of alternatives:choice
## A B C
## 0.48384 0.37505 0.14110
##
## nr method
## 5 iterations, 0h:0m:0s
## g'(-H)^-1g = 2.65E-06
## successive function values within tolerance limits
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## `non-gmo` 1.382253 0.183803 7.5203 5.462e-14 ***
## organic 1.891538 0.187375 10.0950 < 2.2e-16 ***
## salty 0.615630 0.084694 7.2689 3.624e-13 ***
## bitter -0.554799 0.094249 -5.8865 3.945e-09 ***
## nutty -0.686132 0.098393 -6.9734 3.094e-12 ***
## grassy -0.789678 0.102958 -7.6699 1.732e-14 ***
## price -0.047378 0.038500 -1.2306 0.2185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -2254.8As you can see, all of the coefficients are statistically significant except for price. It is likely that how consumers reacted to price is non-linear. The good news is that the price coefficient is negative (remember - people want to pay less, so this coefficient so always be negative).
To calculate WTP, simply take the negative ratio of the parameter for treatment of interest and the price parameter. We can also calculate a confidence interval around the WTP value using the bootstrap or delta method. Those calculations are not presented here.
Random Parameters Logit (RPL)
As we noticed in the conditional logit, price looks to have a heterogeneous effect on edamame choices. We can add consumer heterogeneity directly into the model by allowing the parameters themselves to have a random (distributional) effect. In this case we estimate a mean parameter and a standard deviation of a chosen distribution. Note that the ‘random’ parameter is separate from the random error term of the entire regression. Let’s see what happens if we allow price to have a random component, and will also assume that the parameter is distributed from a normal distribution:
mxl.edamame <- mlogit(choice ~ 0 + `non-gmo` + organic +
salty +bitter + nutty + grassy +
price, panel=FALSE, rpar= c(price = "n") ,
data=edamame.mlogit)
summary(mxl.edamame)##
## Call:
## mlogit(formula = choice ~ 0 + `non-gmo` + organic + salty + bitter +
## nutty + grassy + price, data = edamame.mlogit, rpar = c(price = "n"),
## panel = FALSE)
##
## Frequencies of alternatives:choice
## A B C
## 0.48384 0.37505 0.14110
##
## bfgs method
## 7 iterations, 0h:0m:4s
## g'(-H)^-1g = 3.51E-08
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## `non-gmo` 1.393944 0.194541 7.1653 7.763e-13 ***
## organic 1.903766 0.197044 9.6616 < 2.2e-16 ***
## salty 0.616707 0.086933 7.0940 1.303e-12 ***
## bitter -0.554813 0.090674 -6.1187 9.432e-10 ***
## nutty -0.687273 0.107231 -6.4093 1.462e-10 ***
## grassy -0.790846 0.111132 -7.1162 1.109e-12 ***
## price -0.047074 0.040145 -1.1726 0.241
## sd.price 0.043455 0.079576 0.5461 0.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -2254.7
##
## random coefficients
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## price -Inf -0.07638358 -0.04707357 -0.04707357 -0.01776356 InfThe rpar=() function indicates which parameters are
random and the "n" indicates the normal distribution.
Allowing price to have a random component shows that the mean parameter
value for price is now significant, but there is also significant
variation indicated by the random parameter estimated (sd.price). You
can easily allow the other parameters to be random, and we often do to
allow as much consumer heterogeneity to enter the model as possible.
Mean WTP is calculated as in the conditional logit model. However, if we want to calculate the confidence interval around the WTP value, the calculations become more complicated.
Latent Class Models (LCM)
Another way to allow for heterogeneity is to use a latent class model or clustering methods. The conditional logit model has fallen out of favor in economics literature because it cannot account for heterogeneity, so we have turned to these other methods. Along those lines, we cannot directly include extraneous variables such as gender, age, income, etc. into these models. This is a result of the fact that a single person answers repeated choice questions and the variation across those rows of data has no variation. You can include demographic information by interacting the variables with price or other attributes, but this changes the interpretation of the parameter estimates. As double check your interpretations in the end!
Censored Regression
In this example, we will continue to use edamame data, but from a
different experiment. This was a survey collected by a M.S. student in
School of Plant Sciences as part of their thesis. This survey was done
on paper and the participants are from around Blacksburg and Virginia
Tech campus. In this survey we were interested in how people viewed
alternative end products of edamame (locally grown, organic, fresh, and
on-the-stalk) as compared to the current market option - frozen. The
survey was designed as a one-and-one-half bound contingent valuation
question. This type of experiment is easy to perform on paper surveys
and allows for less biased WTP estimates than a single, dichotomous
choice question. Because the data lends itself to a bounded range for
WTP, we use a censored regression model. This is a bit different as the
dependent variable actually has to parts: a lower and upper bound. The
package used to perform the analysis is the DCchoice
package. Let’s import the data:
edamame_cv_survey <- read_excel("Edamame_CV Example2.xlsx")## New names:
## • `` -> `...4`Depending on how you design your followup questions, depends on the
function used. Since this experiment used a one-and-one-half bound
question format, we use the oohbchoice function.
oohb.edamame <- oohbchoice(formula=Organic_R1 + Organic_R2 ~ 1 + Social_value + Age + Female
| Organic_BL + Organic_BH, data = edamame_cv_survey, na.action = na.omit,
dist="normal")
summary(oohb.edamame) ##
## Call:
## oohbchoice(formula = Organic_R1 + Organic_R2 ~ 1 + Social_value + Age +
## Female | Organic_BL + Organic_BH, data = edamame_cv_survey, na.action = na.omit,
## dist = "normal")
##
## Formula:
## Organic_R1 + Organic_R2 ~ 1 + Social_value + Age + Female | Organic_BL +
## Organic_BH
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.33557 0.63853 -2.0916 0.036471 *
## Social_value 0.22476 0.07738 2.9046 0.003678 **
## Age 0.06463 0.07441 0.8685 0.385101
## Female -0.08887 0.17908 -0.4963 0.619710
## BID 0.22797 0.17260 1.3208 0.186558
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Distribution: normal
## Number of Obs.: 217
## Log-likelihood: -137.932465
##
## LR statistic: 16.976 on 3 DF, p-value: 0.001
## AIC: 285.864929 , BIC: 302.764416
##
## Iterations: 30 8
## Convergence: TRUE
##
## WTP estimates:
## Mean : -4.099098
## Mean : 2.250649 (truncated at the maximum bid)
## Mean : 7.406739 (truncated at the maximum bid with adjustment)
## Median: 2.248344Unlike the conditional and mixed logit models, the parameter estimates are directly interpretative as marginal changes in WTP. This type of analysis is useful as it is easier for non-economist readers to understand from a first glance. However, contingent valuation experiments are continuously criticized for their hypothetical bias and tendencies to overstate WTP estimates.
Conclusions
Throughout this class we have determined sample size, created experimental designs, done data manipulation, and analyzed out data. As I have said throughout, this is just a basic introduction to the world of economic experiments. The methods and survey design can get a lot more complicated, but this should give you the tools to pursue more advanced analysis. This set of notes is also collects information from several packages and sources to host the process all in one place. There is no other place you will find all of this information within one document. I hope that this guide proves useful and guides you on your journey through the world of economic experiments.