Examining fit

This lab will introduce two ways of inspecting the "fit" of a model: posterior mean plots and posterior predictive plots.

Start by loading the data on 2016 annual income, and taking sample of 150 individuals to make our posterior distributions less precise.

d <- read.csv("https://soci620.netlify.com/data/USincome_subset_2016.csv")
# take a random sample of 150 rows
d <- d[sample(nrow(d),150),]
d

A data.frame: 150 × 20

	id	family_size	children_in_hh	female	age	marital_status	times_married	race	hispanic	hispanic_origin	health_ins	health_ins_priv	health_ins_employer	in_school	ed_highschool	ed_college	income	poverty	occ_income	log_income
	<int>	<int>	<int>	<int>	<int>	<chr>	<int>	<chr>	<int>	<chr>	<int>	<int>	<int>	<int>	<int>	<int>	<int>	<lgl>	<int>	<dbl>
606	571365	3	0	0	36	Never married/single	0	White	0	Not Hispanic	0	0	0	0	1	0	25800	FALSE	24	10.158130
2477	2387903	2	0	1	49	Married, spouse present	1	Black/African American/Negro	0	Not Hispanic	1	0	0	0	1	0	13700	FALSE	NA	9.525151
1105	1099431	1	0	0	75	Divorced	2	White	0	Not Hispanic	1	0	0	0	1	0	56600	FALSE	NA	10.943764
88	78464	4	2	0	30	Never married/single	0	White	0	Not Hispanic	1	1	1	0	1	0	30600	FALSE	41	10.328755
1800	1746745	1	0	1	60	Never married/single	0	White	0	Not Hispanic	1	1	1	0	1	0	18100	FALSE	23	9.803667
3134	2990863	4	2	0	43	Married, spouse present	1	Other Asian or Pacific Islander	0	Not Hispanic	1	1	1	0	1	0	85000	FALSE	33	11.350407
2986	2869773	4	2	0	43	Married, spouse present	1	White	0	Not Hispanic	1	1	0	0	1	0	110000	FALSE	26	11.608236
2563	2461937	2	0	0	64	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	21000	FALSE	25	9.952278
3227	3077613	4	2	1	30	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	31000	FALSE	38	10.341742
980	961705	2	0	1	46	Married, spouse present	1	Black/African American/Negro	0	Not Hispanic	1	1	1	0	1	1	125000	FALSE	33	11.736069
2466	2380693	1	0	0	85	Widowed	1	White	0	Not Hispanic	1	0	0	0	1	0	23300	NA	NA	10.056209
2504	2409271	2	0	1	62	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	24000	FALSE	25	10.085809
2289	2210649	1	0	0	60	Divorced	1	White	0	Not Hispanic	1	1	1	0	1	1	38000	FALSE	27	10.545341
3119	2974225	2	0	1	53	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	36000	FALSE	25	10.491274
1911	1856605	2	0	0	36	Married, spouse present	1	White	0	Mexican	1	1	1	0	1	0	30000	FALSE	34	10.308953
1600	1541948	2	0	0	58	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	46800	FALSE	19	10.753638
70	63251	2	0	0	65	Divorced	1	White	0	Not Hispanic	1	0	0	0	1	1	412000	FALSE	62	12.928779
2344	2271390	3	0	1	25	Never married/single	0	American Indian or Alaska Native	0	Not Hispanic	1	0	0	0	1	0	18800	FALSE	13	9.841612
1289	1250200	2	0	0	48	Never married/single	0	White	0	Not Hispanic	1	0	0	0	1	0	13800	FALSE	NA	9.532424
2965	2848590	2	0	1	67	Married, spouse present	2	White	0	Not Hispanic	1	1	1	0	1	0	24400	FALSE	22	10.102338
747	736019	6	3	1	26	Never married/single	0	Black/African American/Negro	0	Not Hispanic	1	0	0	0	1	0	5500	FALSE	22	8.612503
846	828071	4	2	0	43	Never married/single	0	White	0	Not Hispanic	1	1	1	0	1	0	60000	FALSE	33	11.002100
1151	1137032	5	2	1	48	Married, spouse present	3	White	0	Not Hispanic	1	0	0	0	1	0	15800	FALSE	23	9.667765
1013	997966	4	2	1	31	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	95000	FALSE	21	11.461632
435	409170	1	0	0	56	Divorced	1	White	0	Not Hispanic	1	0	0	0	1	0	5000	TRUE	20	8.517193
2870	2775303	4	0	0	49	Never married/single	0	White	0	Mexican	1	0	0	0	1	0	8800	FALSE	NA	9.082507
612	577818	1	0	0	62	Married, spouse absent	1	White	0	Not Hispanic	1	1	1	0	1	1	114500	FALSE	42	11.648330
583	548296	1	0	0	93	Widowed	1	White	0	Not Hispanic	1	1	0	0	1	0	14100	FALSE	NA	9.553930
457	428496	1	0	1	28	Never married/single	0	Two major races	0	Not Hispanic	0	0	0	0	1	0	28800	FALSE	42	10.268131
197	183666	2	0	1	46	Divorced	1	White	0	Not Hispanic	1	1	1	0	1	0	85000	FALSE	33	11.350407
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
2293	2216204	2	0	1	76	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	0	0	7000	FALSE	NA	8.853665
3295	3134474	1	0	1	18	Never married/single	0	White	0	Not Hispanic	1	1	1	1	1	0	2200	NA	13	7.696213
2923	2813298	1	0	1	40	Divorced	1	White	0	Not Hispanic	1	1	1	0	1	1	70000	FALSE	27	11.156251
1284	1247595	5	3	1	33	Married, spouse present	1	White	0	Not Hispanic	1	1	1	1	1	1	46000	FALSE	27	10.736397
3012	2891885	2	0	1	69	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	12000	FALSE	NA	9.392662
2396	2318689	2	0	1	64	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	31800	FALSE	33	10.367222
320	305240	1	0	0	27	Never married/single	0	Other race, nec	0	Mexican	1	1	1	0	1	0	90000	FALSE	42	11.407565
1231	1208900	3	1	1	30	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	0	0	15000	FALSE	22	9.615805
999	982895	6	3	1	46	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	33800	FALSE	43	10.428216
1774	1719562	4	2	1	42	Married, spouse present	1	Japanese	0	Not Hispanic	1	1	1	0	1	0	103000	FALSE	42	11.542484
2667	2577590	1	0	0	49	Never married/single	0	White	0	Not Hispanic	1	1	1	0	1	0	39000	FALSE	38	10.571317
69	62895	2	0	0	51	Married, spouse present	2	White	0	Not Hispanic	0	0	0	0	1	0	30000	FALSE	25	10.308953
1000	984363	4	2	0	48	Married, spouse present	1	Other Asian or Pacific Islander	0	Not Hispanic	1	1	1	0	1	0	132000	FALSE	33	11.790557
446	421992	1	0	1	62	Divorced	1	Black/African American/Negro	0	Not Hispanic	1	1	1	0	1	0	40000	FALSE	16	10.596635
2742	2649201	4	2	0	44	Married, spouse present	1	Black/African American/Negro	0	Not Hispanic	1	1	1	0	1	0	13000	FALSE	42	9.472705
2897	2791523	6	1	1	55	Married, spouse present	1	White	0	Not Hispanic	1	0	0	0	1	0	7000	FALSE	NA	8.853665
1167	1155428	1	0	0	89	Married, spouse absent	1	White	0	Not Hispanic	1	1	1	0	1	1	34000	NA	NA	10.434116
2353	2279738	1	0	0	38	Divorced	1	White	0	Not Hispanic	0	0	0	0	1	0	26000	FALSE	23	10.165852
3009	2886864	1	0	0	27	Never married/single	0	Other Asian or Pacific Islander	0	Not Hispanic	1	1	1	0	1	1	22000	FALSE	42	9.998798
1940	1882346	1	0	0	33	Never married/single	0	Black/African American/Negro	0	Not Hispanic	0	0	0	0	1	0	28000	FALSE	22	10.239960
1586	1529038	5	3	1	34	Married, spouse present	1	White	0	Not Hispanic	1	1	0	0	1	0	24000	FALSE	36	10.085809
2197	2115506	1	0	1	21	Never married/single	0	White	0	Not Hispanic	1	1	1	1	1	0	9000	TRUE	24	9.104980
2311	2231587	4	2	1	54	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	22000	FALSE	25	9.998798
2144	2071582	2	1	0	55	Divorced	2	White	0	Not Hispanic	1	0	0	0	1	0	25000	FALSE	NA	10.126631
19	19751	3	1	0	33	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	1	75000	FALSE	42	11.225243
395	379665	2	1	1	88	Widowed	1	White	0	Not Hispanic	1	1	0	0	1	0	18800	FALSE	NA	9.841612
3276	3118806	3	1	0	56	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	22900	FALSE	19	10.038892
736	722825	4	2	1	46	Married, spouse present	1	Other Asian or Pacific Islander	0	Not Hispanic	1	1	1	0	1	0	33000	FALSE	24	10.404263
2714	2621194	3	1	0	32	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	0	22000	FALSE	16	9.998798
1227	1207810	2	0	0	70	Married, spouse present	1	White	0	Not Hispanic	1	1	1	0	1	1	43800	FALSE	80	10.687389

We can start by building a simple model that predicts a person's log-income using just their age.

library(rethinking)
m1 <- alist(
    log_income ~ dnorm(mu,sigma),
    mu <- a + b1*age,
    a ~ dnorm(8,3),
    b1 ~ dnorm(0,2),
    sigma ~ dunif(0,5)
)
fit1 <- quap(m1,data=d)
# look at the estimates
precis(fit1)

Loading required package: rstan

Loading required package: StanHeaders

Loading required package: ggplot2

rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)

Loading required package: parallel

rethinking (Version 2.13)


Attaching package: ‘rethinking’


The following object is masked from ‘package:stats’:

    rstudent

A precis: 3 × 4

	mean	sd	5.5%	94.5%
	<dbl>	<dbl>	<dbl>	<dbl>
a	9.895820675	0.258221894	9.483132215	10.30850914
b1	0.006979848	0.005069483	-0.001122166	0.01508186
sigma	1.132539930	0.065388349	1.028036719	1.23704314

If the estimate for age seems small, that is because we are measuring age in years and the range of ages in the data is large.

Let's re-estimate using standardized age.

# just to look at the age range in the data
range(d$age)

# construct standardized age
d$age.s <- scale(d$age)
range(d$age.s)

1. 17
2. 93

1. -1.67772886540278
2. 2.48827180233664

# recreate fit1, using standardized age
fit1 <- quap(alist(
    log_income ~ dnorm(mu,sigma),
    mu <- a + b1*age.s,
    a ~ dnorm(8,3),
    b1 ~ dnorm(0,2),
    sigma ~ dunif(0,5)
),data=d)
# look at the estimates
precis(fit1)
exp(fit1@coef[2])

A precis: 3 × 4

	mean	sd	5.5%	94.5%
	<dbl>	<dbl>	<dbl>	<dbl>
a	10.2279588	0.09247971	10.08015832	10.3757592
b1	0.1222032	0.09273354	-0.02600287	0.2704093
sigma	1.1331757	0.06547991	1.02852616	1.2378252

b1: 1.12998372946662

We can start to look at the predictions of our model by plotting our "point estimate" of $\mu$ as a function of standardized age.

# set up plot window for projector visibility
options(repr.plot.width=12, repr.plot.height=8)

# first plot the data points
plot(d$age.s, d$log_income)
# note: here we used the plot(x,y) form of the plot function,
# which plots x on the horizontal axis against y on the vertical.
# There is another way to make the same scatter plot using the
# 'formula' form of plot:
#    plot(log_income~age.s,data=d)
# This takes the form of plot(y~x, data=d). The two forms create
# identical output.

# get the coefficients
cf1 <- coef(fit1)
cf1
a_fit <- cf1['a']
b_fit <- cf1['b1']

# create a grid of hypothetical ages to predict
n_grid <- 400
grid <- seq(min(d$age.s),max(d$age.s),length.out=n_grid)

# plot the "deterministic" portion of the model:
#    mu <- a + b1*age.s

lines(grid , a_fit + b_fit*grid)

a

10.2279587622836

b1

0.122203233918481

sigma

1.13317569469974

This shows us the overall trend predicted by the data, but ignores the uncertainty in our estimates of $\mu$ that comes from uncertainty in a and b1.

One simple way to try to express this uncertainty is to draw the same line as above multiple times, using draws from our posterior rather than the maximum a-posteriori estimates.

# get a sample of coefficients
n_samp <- 300
s <- extract.samples(fit1,n_samp)
head(s)

# plot the data
plot(d$age.s,d$log_income)

# create a grid of hypothetical ages to predict
n_grid <- 400
grid <- seq(min(d$age.s),max(d$age.s),length.out=n_grid)

# create a loop to plot each of the lines
for(i in 1:n_samp){
    # each time we go through this loop, i will take
    # on one of the values between 1 and n_samp (30).
    
    # get the coefficients from the i'th row of the sample s
    a_fit <- s$a[i]
    b_fit <- s$b1[i]
    # plot the line, but make it slightly "translucent" for clarity
    lines(grid , a_fit + b_fit*grid , col="#00000011")
}

A data.frame: 6 × 3

	a	b1	sigma
	<dbl>	<dbl>	<dbl>
1	10.32448	0.17737059	1.1066351
2	10.29114	0.03338418	1.0486046
3	10.26172	0.11108619	1.1285676
4	10.19481	0.02294440	1.0948543
5	10.27887	0.15726755	0.9869483
6	10.13689	0.34857763	1.1873645

This gives us an idea of the spread of uncertainty in our model, but can be hard to read or interpret. A common approach here is to replace the pile of spaghetti in this plot with a shaded region that will show, say, the 90% credible interval for our posterior.

The procedure for this is straightforward:

1. For every point along grid, draw a large number of posterior samples (say 1000) for a and b1.
2. Calculate the value of $\mu$ for each of those samples.
3. Calculate a credible interval on $\mu$ for each point long grid.

There are lots of ways to do this in R, and the rethinking package has some tools to make it easier. But we will first do it "manually" to make it clear what is going on.

# create a grid of hypothetical ages to predict
n_grid <- 400
grid <- seq(min(d$age.s),max(d$age.s),length.out=n_grid)

# how many samples at each grid point?
n_samp <- 1000
s1 <- extract.samples(fit1,n_samp)
head(s1)
# use the `sapply()` function to calculate the 90% credible
# interval at each point x in the grid
mu_ci <- sapply(grid,function(x){
    # this is like a loop, where x will take on every value
    # in the grid
    quantile(s1$a + s1$b1 * x , probs=c(.05,.95))
})
# see what the output of this looks like
mu_ci

# plot the data
plot(d$age.s,d$log_income,pch=16)

# use the `polygon` function to plot a shaded region
polygon(
    x = c(grid,rev(grid)), # rev() returns its argument in reverse order
    y = c(mu_ci[1,],rev(mu_ci[2,])),
    border = NA , col = '#00000055')

A data.frame: 6 × 3

	a	b1	sigma
	<dbl>	<dbl>	<dbl>
1	10.31780	0.04336469	1.223364
2	10.28770	0.05917803	1.002291
3	10.34150	0.15617497	1.207649
4	10.22377	0.30798533	1.060152
5	10.24130	-0.07191926	1.036488
6	10.27759	0.15282044	1.147524

A matrix: 2 × 400 of type dbl

5%	9.741618	9.744198	9.745957	9.748336	9.75127	9.754089	9.756498	9.758661	9.760387	9.761957	⋯	10.13678	10.13619	10.13560	10.13475	10.13386	10.13297	10.13207	10.13118	10.13029	10.12940
95%	10.296917	10.297433	10.297591	10.297695	10.29797	10.298574	10.298795	10.298737	10.298675	10.298585	⋯	10.89463	10.89696	10.89928	10.90159	10.90390	10.90621	10.90855	10.91125	10.91395	10.91664

Looks good! But a bit of a pain. Let's try the same thing using some functions from the rethinking package:

• link() will draw samples from posterior distribution of $\mu$ and calculated predicted values from the data data.
• PI() calculates the percentile interval on a vector.
• shade() is just a convenient wrapper around polygon() that we used above.

# create a grid of hypothetical ages to predict
n_grid <- 400
# create a data frame listing our grid as possible values of age.s
pred_dat <- data.frame(
    age.s = seq(min(d$age.s),max(d$age.s),length.out=n_grid)
)
head(pred_dat)

# how many samples at each grid point?
n_samp <- 1000

# use the 'link' function to take n_samp posterior samples at each point in pred_dat
mu_post <- link(fit1, data=pred_dat, n=n_samp)
# mu_post has one row for each of our posterior samples, and one column
# for each point in our grid.
dim(mu_post)

A data.frame: 6 × 1

	age.s
	<dbl>
1	-1.677729
2	-1.667288
3	-1.656847
4	-1.646406
5	-1.635964
6	-1.625523

1. 1000
2. 400

# calculate the credible interval for every column of mu_post
# `apply()` here will call PI(mu[,i],prob=0.9) for every column i.
mu_ci <- apply(mu_post, 2, PI, prob=0.9)
mu_ci

A matrix: 2 × 400 of type dbl

5%	9.725457	9.728678	9.731687	9.734335	9.736982	9.739629	9.742277	9.744924	9.747572	9.750231	⋯	10.14056	10.14083	10.1411	10.14136	10.14163	10.14189	10.14216	10.14242	10.14269	10.14246
95%	10.307084	10.307032	10.307035	10.307017	10.307045	10.307072	10.307054	10.307005	10.306957	10.306909	⋯	10.89348	10.89654	10.8996	10.90217	10.90457	10.90695	10.90915	10.91148	10.91415	10.91657

# plot the data
plot(d$age.s,d$log_income,pch=16)

# use the `shade()` function from rethinking to draw the shaded region
shade(mu_ci, pred_dat$age.s)