I didn’t mean() to ignore the median()
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
This week’s post follows directly from last week’s
investigation
of data from the 2016 US Census Bureau’s American Community Survey (ACS) Public
Use Microdata Sample
(PUMS). We explored
mean
differences in income across several different types of employment status
(selfemployed, private sector, government, etc.). We found, using bayesian
methods, strong evidence for differences across the categories and were able to
plot them in a variety of formats using ggplot2
and ggsignif
.
While I was very happy with the plot results (which was my primary interest)
something was nagging me the whole time I was writing the post. Our outcome
measure was income
and both our data and the prevailing wisdom from a vast
amount of research is that income data tend to be skewed or suffer from a floor
effect or any number of ways of saying that perhaps the median (or some other
measure of central tendency) is a much better measure than the mean.
Ironically, it is actually very common for income
to be an exemplar in
undergraduate statistics texts to be discussed in context of the dangers of
using the mean. Here’s an example from my favorite undergrad text book Learning
Statistics with R – section
5.1.5.
As you look at our data plot
it is pretty clear that the string of data points on the high end of income in
several categories may be unduly influencing our results. The risk wouldn’t be
quite so bad if the data were consistent across employment categories but a
visual inspection makes us wonder, especially since we can “see” the difference
between the median (the solid black line near the center of the boxplot) and the
mean (plotted here as a dark red box).
For this post then, we’re going to explore comparing medians in a bayesian
framework and along the way spend a minute or two on credible intervals.
Read on! I’ll say it again below but comments and critique are always welcomed
via disqus
or email. You’ll have no trouble finding the icon links in a couple
of places.
To the data cave!
A colleague of mine loves to call her office her “data cave”. I love the term
and it reminds me of the old school Batman series on television, “To the Bat
Cave, Robin”. So let’s go to our virtual data cave and get things setup.
We’ll once again grab the data set from the WinVector Github
site, load the packages we need (hiding the load
messages), set my favorite ggplot theme, and do a little factor cleanup. If you
need a refresher see the last
post
for details.
library(tidyverse) library(BayesFactor) library(ggsignif) library(scales) theme_set(theme_bw()) # set theme location < "https://github.com/WinVector/PDSwR2/raw/master/PUMS/incomedata.rds" incomedata < readRDS(url(location)) incomedata$gp < NULL incomedata$employment < fct_recode(incomedata$employment, "Self Not Inc" = "Self employed not incorporated", "Self Incorporated" = "Self employed incorporated", "Private for Profit" = "Employee of a private for profit", "Private Non Profit" = "Private notforprofit employee", "Federal Government" = "Federal government employee", "State Government" = "State government employee", "Local Government" = "Local government employee" ) incomedata$employment < forcats::fct_reorder( .f = incomedata$employment, .x = incomedata$income, .fun = mean ) # In case I want to reduce down to 3 more fundamental # categories Private, Government, or self employed incomedata$empcategory < fct_collapse(incomedata$employment, "Self" = c("Self Not Inc", "Self Incorporated"), "Private" = c("Private for Profit", "Private Non Profit"), "Government" = c("Federal Government", "State Government", "Local Government") ) str(incomedata)
## 'data.frame': 22241 obs. of 6 variables: ## $ income : num 22000 21000 25800 25000 20200 36000 20000 30000 23000 5000 ... ## $ age : num 24 31 26 27 27 47 24 41 43 21 ... ## $ sex : Factor w/ 2 levels "Male","Female": 1 2 2 2 2 1 1 1 2 2 ... ## $ employment : Factor w/ 7 levels "Self Not Inc",..: 2 4 2 2 4 2 2 2 2 2 ... ## $ education : Factor w/ 9 levels "no high school diploma",..: 1 4 4 6 4 2 1 1 6 2 ... ## $ empcategory: Factor w/ 3 levels "Self","Private",..: 2 2 2 2 2 2 2 2 2 2 ...
Okay we have the data we need. For those who cynically believe that visual
inspection may be an aphorism for “draw any conclusion you like from the
picture” let’s just quickly confirm we have a positive skew with psych::skew
(rounded). In our case the answer is round(psych::skew(incomedata$income), 2)
= 1.78. While we’re at it let’s confirm
my allegation that skewness varies by group using psych::describeBy
. I wasn’t
making it up!
psych::describeBy(incomedata$income, group = incomedata$employment, mat = TRUE, digits = 2) %>% mutate(.data = ., Employment_Type = group1, trimmed_mean = trimmed) %>% select(Employment_Type, mean, sd, median, skew, n, trimmed_mean)
## Employment_Type mean sd median skew n trimmed_mean ## 1 Self Not Inc 41382.50 34209.02 30000 1.73 829 35792.72 ## 2 Private for Profit 51163.99 38744.57 40000 1.82 16170 44932.92 ## 3 State Government 52612.39 29331.94 47150 1.64 1028 49269.63 ## 4 Private Non Profit 52996.20 34389.69 45000 1.85 1573 48085.83 ## 5 Local Government 53720.47 29438.20 48000 1.29 1492 50763.72 ## 6 Federal Government 65311.26 34790.20 60000 0.73 580 62641.60 ## 7 Self Incorporated 66906.44 51983.28 50000 1.52 569 58579.12
We’re skewed! How do we fix it?
There are any number of ways we can address the skew issue (I won’t try methods such
as: hoping your reader/reviewer is inattentive or praying to your favorite
deity). When I say any number I mean a lot of different ways we can address
the issue. While nonexhaustive they generally lump into one of three
categories. We can trim, we can transform, or we can use nonparametric methods
(which is a bit disingenuous since some of the nonparametric methods themselves
rely on transforming our data – we’ll get to that in a minute). Let’s briefly
discuss.

Trim. We can “trim” the high and lows hoping this will eliminate any
troubling values at both ends. You may be surprised to know that capability
is built intomean()
in base R, andmean(incomedata$income, , trim = .1)
will work. For those of you paying attention I grabbed that info when I ran
psych::describeBy
. As you would expect for each of our groups it has
produced a trimmed mean which is closer to the median than our original mean. My
problem with it is threefold. I hate throwing away data. In this case we have
plenty but still. Second, there is no one right answer as to how much we
should trim? Why .1? Why not .05 or .2? It can make you suspicious that we
are simply picking an amount of trimming that best serves our agenda. Finally,
it implicitly makes the extreme values wrong or invalid or “not right”. In
this case there is little evidence that this is true. The income values don’t
look impossible or wrong, they just seem to reflect that there are
people who make much more and that zero is as low as it goes. Personally I
tend to only use trimming when I think the measurement is producing bad
values typically on both ends of the spectrum. 
Transform. One of the common solutions to the income data problem is to
transform the variable. The transformation can take the form of a square
root, a log, or a more complex method such as BoxCox or YeoJohnson. The
later is especially useful if your variable can validly take on negative
values (negative income anyone?). All of these methods are likely to “fix” our
positive skew problem since they will all serve to “pull” our data to the
left in a consistent way. All of them make it more difficult to describe our
results since for many readers the log of income is not a natural concept.
The other limitation here is that there is no one approved transform all will
serve to “improve” our skew but we’re subject to being questioned as to why
we selected one unless it is a disciplinary custom or for reasons of
comparison. 
Nonparametric. For our particular case there are two non parametric methods
we would select. For the overall comparison of income across the 7 levels of
Employment it would be the KruskalWallis. For pairwise comparisons as we did
last post it would be the MannWhitney U (a.k.a. Wilcoxon rank sum) test. If
you look back you’ll even notice that the Wilcoxon rank sum test
wilcox.test()
is the default forggsignif
. As I mentioned earlier these
tests are actually a sort of transformation as well since they rely on the
ranks of our data rather than the data itself. I admit to having a bias
towards (maybe even a fondness for) nonparametric tests. While they aren’t
perfect and have limitations I take a certain amount of solace in knowing
that any supported findings were by using tests that at least eliminate
some number of assumptions.
So after last week’s posting I went sleuthing for the Bayesian equivalent to the
Mann Whitney and KruskalWallis tests. Yes I recognize that with powerful
software like Stan you can pretty much construct any data analysis you want but
I was hoping for something more “off the shelf” something very much like
BayesFactor::ttestBF
. No luck, I won’t claim I searched everywhere but I did
search a lot for anything already done for R. Lot’s of theoretical papers and
some questions but no solutions. If I missed it, and someone knows about it,
please share.
Since I have been tracking the progress of JASP
(unsolicited praise – it
looks like a great piece of software especially if you come from the world of
SPSS
which I do). A lot of the calculation that is done under the very nice
GUI is powered by R code and libraries. I noticed that they had recently
implemented a Bayesian version of the Mann Whitney so I decided to give it a
try. No joy! Using the very same data I have been using for these blog posts
it ran for well over an hour and hadn’t completed. When I came back the next day
it had completed and produced what appeared to be reasonable results on the
surface but I had to do some sleuthing because the very long run time had me
confused.
The underlying code is open source and you can see it
here.
After a little bit of fiddling I was able to separate the JASP interface pieces
from the R code that produces the results. I was thwarted again. Even on a much
smaller dataset like ggplot2::mpg the results are slow and unstable. That’s not
meant as a complaint I know they are working on it and I know working through
the ties is challenging.
A Eureka and a “YouRankA”
Temporarily thwarted – something in the back of my mind made me remember dimly
from the past considering the case of t tests on rankings. It had been a long
time (nowadays it so trivial to run for example wilcox.test
) and I wasn’t quite
sure how comparable the results would be so I went looking. Found research on
the frequentist side at least that the methodologies produce consistent results
(Zimmerman, 2012).
So if I can’t have a mannwhitneyBF
maybe I can ttestBF(rank(income))
? Seemed
promising and the initial results seemed to be reasonable as well. The only
challenge seemed to be that the distribution of rank(anything)
tends to be
anything but normal. The result tends to be very very flat with almost no tails
at all. Almost a uniform distribution. Hmmmmm, which is a worse violation of
ttest assumptions, skewed or flat?
As I was puzzling on that I stumbled upon another avenue of exploration from a
blog post by Andrew Gelman at Columbia
here
and
here.
Quoting…
“Takehome message: I am not saying that the ranktheninversenormaltransform
strategy is always or even often a good idea. What I’m saying is that, if you
were planning to do a rank transformation before analyzing your data, I
recommend this zscore approach rather than the classical Wilcoxon method.” Gelman (July 2015)
That sounded promising. Was there an easy way for someone like me to implement
his better than Wilcoxon suggestion (point taken that he repeatedly says “I’m
not saying you should transform the data. I’m saying that if you were going to
do Wilcoxon, then you’re already doing a rank transformation, in which case I
think it makes sense to do it in an expandable way, rather than getting hung up
on some pretabulated test statistics.”)?
Turns out there is. And it’s easy. And it’s already written and packaged in the
RNOmni
package on CRAN. The
function is rankNorm
and straight from the description what it does is:
“Applies the rank based inverse normal transform (INT) to a numeric vector. The
INT can be broken down into a twostep procedure. In the first, the
observations are transformed onto the probability scale using the empirical cumulative
distribution function (ECDF). In the second, the observations are transformed
onto the real line, as Zscores, using the probit function.”. The key component
of the function is qnorm((r  k)/(n  2 * k + 1))
where r is the rank, n the
number of observations, and k = 3/8 = the Blom transform.
I just realized I’ve been rambling for quite some time without a single graphic
or plot. Perhaps this is a good time to just display the original income data
and what the various transformations would do to them. I’ll plot them as density
plots and superimpose the theoretical normal curve in a different color.
Let’s use cowplot
and make 7 plots in order:
 The original
income
data  Transform income to a z score (mean = 0, sd = 1)
scale(income)
 Transform
sqrt(income)
 Transform
log(income)
 Transform
VGAM::yeo.johnson(income, lambda = .1)
 Transform
rank(income)
 Transform
RNOmni::rankNorm(income)
So as a first step we’ll use the mutate
function to create these new variables. We’ll deliberately name them all so they end with “_income”. In a minute we’ll use that little bit of consistency to make our lives easier.
incomedata < incomedata %>% mutate(z_income = as.numeric(scale(income)), sqrt_income = sqrt(income), log_income = log(income), yj_income = VGAM::yeo.johnson(income, lambda = .1), rank_income = rank(income), rint_income = RNOmni::rankNorm(income) ) glimpse(incomedata)
## Observations: 22,241 ## Variables: 12 ## $ income <dbl> 22000, 21000, 25800, 25000, 20200, 36000, 20000, 300… ## $ age <dbl> 24, 31, 26, 27, 27, 47, 24, 41, 43, 21, 29, 30, 28, … ## $ sex <fct> Male, Female, Female, Female, Female, Male, Male, Ma… ## $ employment <fct> Private for Profit, Private Non Profit, Private for … ## $ education <fct> "no high school diploma", "some college credit, no d… ## $ empcategory <fct> Private, Private, Private, Private, Private, Private… ## $ z_income <dbl> 0.791115497, 0.817539654, 0.690703699, 0.7118430… ## $ sqrt_income <dbl> 148.32397, 144.91377, 160.62378, 158.11388, 142.1267… ## $ log_income <dbl> 9.998798, 9.952278, 10.158130, 10.126631, 9.913438, … ## $ yj_income <dbl> 17.17967, 17.05353, 17.61618, 17.52934, 16.94867, 18… ## $ rank_income <dbl> 3941.5, 3694.5, 5444.5, 5082.0, 3531.0, 9310.0, 3191… ## $ rint_income <dbl> 0.9260927, 0.9697185, 0.6910214, 0.7438702, 0.9…
A quick glimpse
indicates we have what we need. So we have seven variables and
we’re going to create 7 identical plots using ggplot
. For each plot we’ll
plot the density
function for our observations (a more granular histogram so
to speak). We’ll also plot the theoretical normal curve for that variable using
dnorm()
with the mean and standard deviation for the variable.
We could just write (mainly cut and paste) the call to ggplot
7 times. But
that seems wasteful. Instead we’ll use purrr:pmap
to do what we want with a
lot less code. pmap
wants a list that has as it’s components the 3 things that
are going to vary for each of the seven iterations. The 3 things are the name of
the variable, the mean, and the sd for the variable. Since ggplot
can be
fickle about passing it bare variable names we’ll pass it character strings for
the name. The mean and the sd of course are numeric types. Again we could
manually build those vectors but that would be tedious so let’s use some
dplyr::select_at
and summarise_at
statements to get what we need. Since we
were careful to name our variables so that they all end with “income” we can use
vars(ends_with("income"))
to grab just what we need and not all the variables
in incomedata
.
a < incomedata %>% select_at(vars(ends_with("income"))) %>% names b < incomedata %>% summarise_at(vars(ends_with("income")), mean, na.rm = TRUE) c < incomedata %>% summarise_at(vars(ends_with("income")), sd, na.rm = TRUE) plot_list < list( which = a, means = b, stddevs = c ) plot_list
## $which ## [1] "income" "z_income" "sqrt_income" "log_income" "yj_income" ## [6] "rank_income" "rint_income" ## ## $means ## income z_income sqrt_income log_income yj_income rank_income ## 1 51939.1 2.884567e17 214.8975 10.61248 18.97702 11121 ## rint_income ## 1 4.205901e05 ## ## $stddevs ## income z_income sqrt_income log_income yj_income rank_income ## 1 37844.16 1 75.88434 0.7348988 2.09446 6419.503 ## rint_income ## 1 0.9991945
That took more text to describe than it took code to enact!
Now we invoke ggplot
as a function .f
. Inside the function call you’ll see ..1, ..2 & ..3 those mark the places where pmap
will substitute in the appropriate value of plot_list
. The output of pmap
is itself a list. A list of ggplot
objects. We’ll take that list and give it to cowplot::plot_grid
along with some labeling and formatting and voila!
list_of_plots < pmap(.l = plot_list, .f = ~ ggplot(data = incomedata, aes_string(x = ..1) ) + geom_density(alpha = .2) + stat_function(fun = dnorm, color = "red", args = list(mean = ..2, sd = ..3) ) ) cowplot::plot_grid(plotlist = list_of_plots, labels = c("Income (Original)", "Z score Income", "SQRT Income", "Log Income", "Yeo Johnson Income", "Ranked Income", "Rank Inv Norm Income"), label_size = 12, vjust = 1, hjust = 0, scale = .9)
Not surprisingly the rank inverse normal transform works as advertised. We have
succeeded in taking the ranks and backing them into a near perfect normal curve.
That’s good but we did that work to prepare for comparing across employment
groups so we would be wise to take a look at our transformed income variable
rint_income
by group. Once again we’ll make use of psych::describeBy
. This
time we’ll also include the standard deviation since we’d like to have some
sense of whether the variances are at least roughly equal among the groups. While
we’re at it we’ll plot the density curves on a by group basis.
psych::describeBy(incomedata$rint_income, group = incomedata$employment, mat = TRUE, digits = 2) %>% mutate(.data = ., Employment_Type = group1, trimmed_mean = trimmed) %>% select(Employment_Type, mean, sd, median, skew, n, trimmed_mean)
## Employment_Type mean sd median skew n trimmed_mean ## 1 Self Not Inc 0.40 1.10 0.49 0.04 829 0.41 ## 2 Private for Profit 0.04 1.02 0.05 0.09 16170 0.05 ## 3 State Government 0.14 0.80 0.18 0.29 1028 0.17 ## 4 Private Non Profit 0.10 0.89 0.12 0.06 1573 0.10 ## 5 Local Government 0.16 0.83 0.20 0.51 1492 0.20 ## 6 Federal Government 0.43 0.88 0.54 0.71 580 0.49 ## 7 Self Incorporated 0.31 1.13 0.27 0.02 569 0.31
incomedata %>% ggplot(aes(x = rint_income, fill = employment, color = employment) ) + geom_density(alpha = .2)
Well, we certainly seem to have improved the distributions and made them more
suitable for parametric manipulation. Certainly looks like there are shifts in
central tendency by group. For the astute observer yes I’ll acknowledge that the
variances among groups aren’t “equal”. But an F_{max} ratio of
2 isn’t that awful.
From last post let’s bring back our comparisons_list
function which will
create a list of all the possible paired comparisons for us. We’ll store them in a
list called comp.list
.
comparisons_list < function(data, x) { data < dplyr::select( .data = data, x = !!rlang::enquo(x) ) %>% dplyr::mutate(.data = ., x = droplevels(as.factor(x))) grouplevels < levels(data$x) g1_list < combn(grouplevels, 2)[1, ] g2_list < combn(grouplevels, 2)[2, ] comparisons_list < lapply( 1:length(g1_list), function(i) c( combn(grouplevels, 2)[2, i], combn(grouplevels, 2)[1, i] ) ) return(comparisons_list) } comp.list < comparisons_list(incomedata, employment)
Welcome back my friends to the show that never ends (ELP)
Last week we produced this plot using the original income data and using
ggsignif
. We ran the wilcox.test
against all the pairings to produce p
values on the plot. Here it is unretouched to refresh your memory.
If Gelman and Zimmerman have steered me correctly (or lol if I paid enough
attention) we should be able to substitute in rint_income
our ranked inverse
normalized transformed version of income
and get very similar results using
the very same code and wilcox.test
.
ggplot(data = incomedata, aes( x = employment, y = rint_income, fill = employment, group = employment )) + geom_boxplot(show.legend = FALSE) + geom_signif( comparisons = comp.list, step_increase = .1 ) + scale_y_continuous(breaks = seq(from = 4, to = 4, by = 1) ) + ggtitle(label = "ACS 2016 Rank Inverse Transformed Income by Employer Type", subtitle = "Mann Whitney multiple comparisons non directional hypothesis using wilcox.test")
ggplot(data = incomedata, aes( x = employment, y = rank_income, fill = employment, group = employment )) + geom_boxplot(show.legend = FALSE) + geom_signif( comparisons = comp.list, step_increase = .1 ) + scale_y_continuous(breaks = seq(from = 0, to = 25000, by = 5000) ) + ggtitle(label = "ACS 2016 Ranked Income by Employer Type", subtitle = "Mann Whitney multiple comparisons non directional hypothesis using wilcox.test")
Reminder – Note that in the twosample case the estimator for the difference in location parameters does not estimate the difference in medians (a common misconception) but rather the median of the difference between a sample from x and a sample from y.
Very similar indeed! Notice that the reported p values are essentially identical
even as we can see from the box plots that the transformations are having an
effect.
But of course this has all been prep work to this point. What I’m really after
is to shift over to a bayesian framework and generate bayes factors that are the
equivalent of the frequentist’s Mann Whitney. So let’s get on with that!
Not having been able to find the R code to do the work directly I’ve transformed
the income variable and will now apply the ttestBF
function. Last
post
I wrote a little function called pairwise_bf
that generates the bayes factors
and manipulates them in a way that makes them suitable for plotting with ggplot
and ggsignif
. Here it is again.
pairwise_bf < function(x = NULL, y = NULL, display_type = "bf", k = 2) { results < ttestBF(x = x, y = y) %>% extractBF() %>% mutate(support = case_when( bf < .01 ~ "extreme BF01", bf < .03 & bf >= .01 ~ "very strong BF01", bf < .1 & bf >= .03 ~ "strong BF01", bf < 1 / 3 & bf >= .1 ~ "moderate BF01", bf < 1 & bf >= 1 / 3 ~ "anecdotal BF01", bf >= 1 & bf < 3 ~ "anecdotal BF10", bf >= 3 & bf < 10 ~ "moderate BF10", bf >= 10 & bf < 30 ~ "strong BF10", bf >= 30 & bf < 100 ~ "very strong BF10", bf >= 100 ~ "extreme BF10" )) %>% mutate(logged = case_when( bf < 1 ~ paste("log(BF01) = ", round(log(1 / bf), k)), bf >= 1 ~ paste("log(BF10) = ", round(log(bf), k)) )) %>% mutate(human = case_when( bf < .000001 ~ "BF01 >= 1,000,000 : 1", bf < .001 & bf >= .000001 ~ "BF01 >= 1,000 : 1", bf < .01 & bf >= .001 ~ "BF01 >= 100 : 1", bf < 1 & bf >= .01 ~ paste("BF01 = ", round(1 / bf, k), ": 1"), bf >= 1 & bf < 100 ~ paste("BF01 = ", round(bf, k), ": 1"), bf >= 100 & bf < 1000 ~ "BF10 >= 100 : 1", bf >= 1000 & bf < 1000000 ~ "BF10 >= 1,000 : 1", bf >= 1000000 ~ "BF10 >= 1,000,000 : 1" )) if (display_type == "support") { results < mutate(results, p.value = support) } else if (display_type == "log") { results < mutate(results, p.value = logged) } else if (display_type == "human") { results < mutate(results, p.value = human) } else { results < mutate(results, p.value = bf) } return(results) } # pairwise_bf(incomedata$employment, incomedata$rint_income) comp.list2 < comparisons_list(incomedata, empcategory) theme_set(hrbrthemes::theme_ipsum())
Using the employment
factor in the original dataset has been useful so far.
Having 7 levels has allowed us to make sure that the functions supported lots of
comparisons. The downside is the resulting plots have an overwhelming amount of
information. So let’s shift over to the variable empcategory
which collapses
employment
to just 3 levels. We’ll be comparing the selfemployed, to those
employed in the private sector, to those employed in local, state or federal
jobs.
Seems reasonable to believe that the income levels across those broad
categories, might be different. I don’t have any strong prior beliefs or
information before this data. It’s not an area I study or have a lot of prior
information amount. Clearly, it’s out there, the ACS survey has been running for
years, but for now I’m perfectly happy to say my priors are flat.
What we’re going to do:

Show all the possible pairwise comparisons for
empcategory
as bayes factors inggplot
(we’ve already built the list above incomp.list2
) 
Change the theme just to make the plots nicer
theme_set(hrbrthemes::theme_ipsum())

Apply the
ttestBF
test from theBayesFactor
package to the data piped
intoggplot
. We’ll compareincome
(the original untransformed data),
rint_income
which is our rank inverse transformed variant, and
rank_income
which is the income data with a simple rank transformation.
N.B. – I’ve expressed the ranks so they are directionally equivalent to the
raw data so the highest recorded income in our data is $250,000 which
corresponds to rank of 22226.5 and after inverse transformation = 3.20299 
Display the BF_{10} and/or BF_{01} values that are between 1 and 100 as
is (rounded to 2 digits) and then create ranges between 100 to 1,000, 1,000 to one million and
greater than one million. At some point talking about odds over 100:1 (in my
humble opinion) loses the need for precision. After all is there really much
difference between odds of 1,000:1 versus 1,001:1? 
Since this post is very much about
mean()
versusmedian()
both of them
are plotted on top of the violin geom. The mean is displayed as a black box
and the median is displayed as a dark red circle.
Caveat Emptor
Before I draw any conclusions about the results a reminder. While the data are
“real” drawn from the U.S. Census Bureau they most assuredly are not totally
representative of the entire U.S. population! The very fact that the max value
for income is $250,000 tells us that this isn’t everyone and I’m not sure at all
how “representative” it is. Please don’t stray from my desire to investigate the
statistical properties of the methods to think I’m trying to convince you of
some conclusion about the income levels of Americans.
Conclusions we can draw.

Distributions matter! There is a difference between investigating the mean()
and the median() income and using bayesian methods doesn’t change that. 
Bayes factors give us information in both directions! When you look at the
first plot using the mean for untransformed income values BF_{01} = 29.36
tells you that the evidence is “strong” that the mean income for the
selfemployed is the same as those in the private sector. 
Whether we use simple rank transformation or reverse inverse transformed
ranks our overall conclusions are similar. We have weaker “anecdotal”
evidence that there is no difference between the selfemployed and those in
the private sector. To say it differently we are less convinced from our data
that there is no difference. On the other hand we have beyond very strong
evidence that government employee income is different than either the self
employed or those in the private sector. 
The difference in the bayes factors generated by using ranks versus rank
inverse transformation for our data were not substantively important in
magnitude. Yes they were different, for example 1.63 versus 2.48, but don’t
lead us to any practically different conclusions. 
Not withstanding #4 above until there is a “true” equivalent to the Mann
Whitney for bayesian analysis I personally am far more comfortable applying a
“ttest” like analysis to the data distribution that results from RINT.
Let’s get credible
I know this post is already pretty long but I promised myself I’d capture what I
learned about generating credible intervals as part of this workflow. I had been
asked the question by a reader I did the work and posted it in a comment but now
is my chance to capture it more formally.
In my usual fashion I’m not going to try and explain the concept of a credible
interval in general when others have done so already and probably done a better
job. So if you’re unfamiliar with the concept I recommend a trip to
Wikipedia to read up. A simple
Google search will also yield lots of informative answers.
What I will focus on for the remainder of this post is simply helping you with
the mechanics of constructing credible intervals if you have been following this
series of posts.
So far we have been making heavy use of BayesFactor::ttestBF
to accomplish our
objectives. For example continuing with our last example our plot says the
difference between the mean income for the privately employed and government
workers is “BF10 >= 1,000:1”. We got that nice neat answer as a result of:
library(bayestestR) tempdf < incomedata %>% filter(empcategory == "Private"  empcategory == "Government") %>% droplevels %>% as.data.frame ttestBF(formula = income ~ empcategory, data = tempdf)
## Bayes factor analysis ##  ## [1] Alt., r=0.707 : 357043.1 ±0% ## ## Against denominator: ## Null, mu1mu2 = 0 ##  ## Bayes factor type: BFindepSample, JZS
There it is – 357,043:1 to be more precise. And note the part that says
“Null, mu1mu2 = 0”. That looks familiar that’s a t test alright it is the
difference between the means = 0. So let’s take the output from that ttestBF
and but it in an object called bf.object
. Sadly it’s a rather complex S4
object that doesn’t contain exactly what we need so a little more work is
required. What we need to do is run some Markov chain Monte Carlo
methods to sample from
the posterior distribution to get our answer. While that sounds wickedly complex
the function we need is built right into BayesFactor
as posterior
so one
line does the trick posterior.data < posterior(bf.object, iterations = 10000)
. What’s in posterior.data
? About 10,000 rows so let’s head
it.
The first column mu
is the overall mean, the second column beta (Private  Government)
is exactly what we need to work on it’s the difference in means
“mu1mu2” between Private and Government. The sign tells us that at least
for the first 7 rows Government is higher which agrees with our plot. I have no
desire to scroll through 10,000 rows so let’s get some summary information and
even a plot of our 10,000 samples.
bf.object < ttestBF(formula = income ~ empcategory, data = tempdf) posterior.data < posterior(bf.object, iterations = 10000) head(posterior.data)
## Markov Chain Monte Carlo (MCMC) output: ## Start = 1 ## End = 7 ## Thinning interval = 1 ## mu beta (Private  Government) sig2 delta g ## [1,] 53712.39 4138.397 1388783139 0.11104907 0.3962511 ## [2,] 53328.80 3951.455 1372469541 0.10666102 0.5661250 ## [3,] 53798.28 3568.424 1396930789 0.09547487 0.1124928 ## [4,] 53216.72 3013.836 1401327895 0.08051000 0.2242916 ## [5,] 52885.07 2276.827 1415166246 0.06052381 0.4882304 ## [6,] 52591.16 3624.291 1381320679 0.09751598 0.5180922 ## [7,] 53115.20 3875.964 1410243480 0.10321260 2.2364977
summary(posterior.data)
## ## Iterations = 1:10000 ## Thinning interval = 1 ## Number of chains = 1 ## Sample size per chain = 10000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Timeseries SE ## mu 5.343e+04 3.639e+02 3.639e+00 6.213e+00 ## beta (Private  Government) 4.197e+03 7.303e+02 7.303e+00 1.269e+01 ## sig2 1.396e+09 1.361e+07 1.361e+05 1.361e+05 ## delta 1.124e01 1.955e02 1.955e04 3.396e04 ## g 2.744e+00 4.663e+01 4.663e01 4.663e01 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% ## mu 5.271e+04 5.318e+04 5.343e+04 5.367e+04 ## beta (Private  Government) 5.614e+03 4.689e+03 4.198e+03 3.718e+03 ## sig2 1.370e+09 1.386e+09 1.395e+09 1.405e+09 ## delta 1.501e01 1.255e01 1.123e01 9.948e02 ## g 6.897e02 1.877e01 3.776e01 9.061e01 ## 97.5% ## mu 5.413e+04 ## beta (Private  Government) 2.745e+03 ## sig2 1.422e+09 ## delta 7.344e02 ## g 1.009e+01
plot(posterior.data[,"beta (Private  Government)"])
Whether we look at the right hand plot, or the tabular information Mean =
4.188e+03, 50% = 4.182e+03, sure looks like the difference estimate is about ~
$4,180. Also telling is 2.5% = 5.615e+03 and 97.5% = 2.775e+03 which means
95% of our estimates were between $5,615 and $2,775! Anyone care to guess what
the 95% credible interval is?
Well of course we don’t have to guess. We can take our posterior.data
object and change it to a dataframe. After we’ve done that we can use bayestestR::ci
(you may have noticed I loaded the library a few lines back). And there you have it the 89% credible interval is [5352.79, 3042.40]. Wait. What? 89%????
posterior.data < as.data.frame(posterior.data) str(posterior.data)
## 'data.frame': 10000 obs. of 5 variables: ## $ mu : num 53712 53329 53798 53217 52885 ... ## $ beta (Private  Government): num 4138 3951 3568 3014 2277 ... ## $ sig2 : num 1.39e+09 1.37e+09 1.40e+09 1.40e+09 1.42e+09 ... ## $ delta : num 0.111 0.1067 0.0955 0.0805 0.0605 ... ## $ g : num 0.396 0.566 0.112 0.224 0.488 ...
bayestestR::ci(posterior.data$`beta (Private  Government)`)
## # Credible Interval ## ## 89% CI ## [5349.28, 3016.02]
For an explanation of why please see the bayestestR
documentation.
While you’re there good chance to review what a credible interval is, how it is
interpreted, and our next topic which is that there are at least two equally
valid ways of computing it. The first we used already it’s the Equaltailed
Interval (ETI) which is the default when you invoke ci
. The other is Highest
Density Interval (HDI) hdi()
. They do a good job of
explaining
so I won’t repeat merely show you how.
ci()
will accept a vector with multiple cutoffs so let’s do 89% and 95% once
with equaltailed and then with highest density. To be fair at least for the
equal tailed you can use base R’s quantile
function to get the same answers.
BTW back a little bit when we were looking at summary(posterior.data)
the
answer for 95% was right there just hard to pick out.
ci(posterior.data$`beta (Private  Government)`, ci = c(.89, .95))
## # Credible Intervals ## ## 89% CI ## [5349.28, 3016.02] ## ## ## 95% CI ## [5614.50, 2745.04]
hdi(posterior.data$`beta (Private  Government)`, ci = c(.89, .95))
## # Highest Density Intervals ## ## 89% HDI ## [5327.73, 2999.18] ## ## ## 95% HDI ## [5625.08, 2759.35]
round(quantile(posterior.data$`beta (Private  Government)`, probs = c(.025, .055, (1  .055), (1  .025) ) ), 2)
## 2.5% 5.5% 94.5% 97.5% ## 5614.50 5349.28 3016.02 2745.04
So given our observed data, the mean difference in income between the private sector respondents and the government employed respondents has a 95% probability of falling between $5,614.99 and $2,774.84.
We can preform the exact same operations on rint_income
if we want to convince
ourselves that the mean difference (which for rint_income
is more akin to the
median difference) is not equal to zero). The problem of course is that it
doesn’t return a dollar value and our transformations are nonlinear.
bf.object < ttestBF(formula = rint_income ~ empcategory, data = tempdf) posterior.data < posterior(bf.object, iterations = 10000) summary(posterior.data)
## ## Iterations = 1:10000 ## Thinning interval = 1 ## Number of chains = 1 ## Sample size per chain = 10000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Timeseries SE ## mu 0.08782 0.009480 9.480e05 1.581e04 ## beta (Private  Government) 0.22912 0.019066 1.907e04 3.165e04 ## sig2 0.96652 0.009474 9.474e05 9.474e05 ## delta 0.23306 0.019439 1.944e04 3.222e04 ## g 3.00664 61.759312 6.176e01 6.176e01 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## mu 0.06913 0.08155 0.08785 0.09424 0.1062 ## beta (Private  Government) 0.26596 0.24200 0.22917 0.21659 0.1917 ## sig2 0.94811 0.96013 0.96647 0.97285 0.9851 ## delta 0.27039 0.24625 0.23310 0.22022 0.1948 ## g 0.07617 0.19808 0.40026 0.94421 9.5910
plot(posterior.data[,"beta (Private  Government)"])
posterior.data < as.data.frame(posterior.data) bayestestR::ci(posterior.data$`beta (Private  Government)`)
## # Credible Interval ## ## 89% CI ## [0.26, 0.20]
ci(posterior.data$`beta (Private  Government)`, ci = c(.89, .95))
## # Credible Intervals ## ## 89% CI ## [0.26, 0.20] ## ## ## 95% CI ## [0.27, 0.19]
hdi(posterior.data$`beta (Private  Government)`, ci = c(.89, .95))
## # Highest Density Intervals ## ## 89% HDI ## [0.26, 0.20] ## ## ## 95% HDI ## [0.27, 0.19]
A final look at why it matters
So I wanted to end this post by reinforcing why mean()
vs median()
matters.
Continuing to focus on the differences between income in the private sector and
government employees remember we saw the BF associated with the difference
between the two groups increase an order of magnitude when we shifted from
assessing mean differences to median differences. Why?
Well take a look at the differences expressed in the original dollar units. The
mean difference is about ~$4,200, the median difference is more than twice as
much at ~$10,000, that’s a big difference.
psych::describeBy(incomedata$income, group = incomedata$empcategory, mat = TRUE, digits = 2) %>% mutate(.data = ., Employment_Type = group1, trimmed_mean = trimmed) %>% select(Employment_Type, mean, sd, median, skew, n, trimmed_mean)
## Employment_Type mean sd median skew n trimmed_mean ## 1 Self 51771.00 44154.93 40000 1.82 1398 44342.71 ## 2 Private 51326.42 38381.11 40000 1.82 17743 45221.13 ## 3 Government 55521.62 30830.33 50000 1.28 3100 52263.65
a < median(incomedata[incomedata$empcategory == "Private", "income"]) b < median(incomedata[incomedata$empcategory == "Government", "income"]) c < mean(incomedata[incomedata$empcategory == "Private", "income"]) d < mean(incomedata[incomedata$empcategory == "Government", "income"]) median_diff_priv_govt < a  b median_diff_priv_govt
## [1] 10000
mean_diff_priv_govt < c  d mean_diff_priv_govt
## [1] 4195.196
For a final plot we’ll show it in the original units and then to make it easier
to see give you a zoomed perspective with coord_cartesian()
. The dashed lines
are the median values by group (which are colorcoded) and the solid lines are
the mean.
tempdf %>% ggplot(aes(x = income, fill = empcategory, color = empcategory) ) + geom_density(alpha = .3) + geom_vline(xintercept = a, linetype = 5, color = "coral1") + geom_vline(xintercept = b, linetype = 5, color = "cyan3") + geom_vline(xintercept = c, color = "coral1") + geom_vline(xintercept = d, color = "cyan3") + scale_x_continuous(label = dollar)
tempdf %>% ggplot(aes(x = income, fill = empcategory, color = empcategory) ) + geom_density(alpha = .3) + geom_vline(xintercept = a, linetype = 5, color = "coral1") + geom_vline(xintercept = b, linetype = 5, color = "cyan3") + geom_vline(xintercept = c, color = "coral1") + geom_vline(xintercept = d, color = "cyan3") + scale_x_continuous(label = dollar) + coord_cartesian(xlim = c(10000, 100000))
Done
Hope you enjoyed the post. Comments always welcomed.
Chuck
This work is licensed under a Creative Commons AttributionShareAlike 4.0 International License
Rbloggers.com offers daily email 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/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.