Gaussian multiple linear regression

Prior predictive for model comparison

1 Dataset

The data is a 1976 Panel Study of Income Dynamics, based on data for the previous year, 1975. Of the 753 observations, the first 428 are for women with positive hours worked in 1975, while the remaining 325 observations are for women who did not work for pay in 1975. A more complete discussion of the data is found in Mroz [1987], Appendix 1. Thomas A. Mroz (1987) The Sensitivity of an Empirical Model of Married Women’s Hours of Work to Economic and Statistical Assumptions. Econometrica, Vol. 55, No. 4 (July 1987), pp. 765-799. Stable URL: http://www.jstor.org/stable/1911029. The variables in the dataset are as follows:

• LFP “A dummy variable = 1 if woman worked in 1975, else 0”;

• WHRS “Wife’s hours of work in 1975”;

• KL6 “Number of children less than 6 years old in household”;

• K618 “Number of children between ages 6 and 18 in household”;

• WA “Wife’s age”;

• WE “Wife’s educational attainment, in years”;

• WW “Wife’s average hourly earnings, in 1975 dollars”;

• RPWG “Wife’s wage reported at the time of the 1976 interview (not the same as the 1975 estimated wage). To use the subsample with this wage, one needs to select 1975 workers with LFP=1, then select only those women with non-zero RPWG. Only 325 women work in 1975 and have a non-zero RPWG in 1976.”;

• HHRS “Husband’s hours worked in 1975”;

• HA “Husband’s age”;

• HE “Husband’s educational attainment, in years”;

• HW “Husband’s wage, in 1975 dollars”;

• FAMINC “Family income, in 1975 dollars. This variable is used to construct the property income variable.”;

• MTR “This is the marginal tax rate facing the wife, and is taken from published federal tax tables (state and local income taxes are excluded). The taxable income on which this tax rate is calculated includes Social Security, if applicable to wife.”;

• WMED “Wife’s mother’s educational attainment, in years”;

• WFED “Wife’s father’s educational attainment, in years”;

• UN “Unemployment rate in county of residence, in percentage points. This taken from bracketed ranges.”;

• CIT “Dummy variable = 1 if live in large city (SMSA), else 0”;

• AX “Actual years of wife’s previous labor market experience”;

2 Importing the data

data = read.table("http://hedibert.org/wp-content/uploads/2020/01/mroz-data.txt",header=TRUE)

attach(data)

names = c("LFP","WHRS","KL6","K618","WA","WE","WW","RPWG","HHRS",
          "HA","HE","HW","MTR","WMED","WFED","UN","CIT","AX")

n = nrow(data)
y = scale(log(FAMINC))
X = cbind(LFP,WHRS,KL6,K618,WA,WE,WW,RPWG,HHRS,HA,HE,HW,MTR,WMED,WFED,UN,CIT,AX)
X = scale(X)
p = ncol(X)

3 Maximum likelihood estimation (MLE)

ols = lm(y~X-1)
summary(ols)

## 
## Call:
## lm(formula = y ~ X - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5881 -0.1919  0.0090  0.1870  3.0093 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## XLFP  -0.088521   0.029839  -2.967  0.00311 ** 
## XWHRS  0.132256   0.028473   4.645 4.03e-06 ***
## XKL6   0.019941   0.019610   1.017  0.30954    
## XK618  0.114282   0.018957   6.029 2.62e-09 ***
## XWA    0.055325   0.038119   1.451  0.14711    
## XWE   -0.005337   0.023802  -0.224  0.82265    
## XWW    0.062273   0.025077   2.483  0.01324 *  
## XRPWG  0.036699   0.025360   1.447  0.14829    
## XHHRS  0.059383   0.021027   2.824  0.00487 ** 
## XHA    0.041883   0.036653   1.143  0.25354    
## XHE    0.038257   0.022749   1.682  0.09306 .  
## XHW    0.147563   0.031960   4.617 4.59e-06 ***
## XMTR  -0.711872   0.032891 -21.643  < 2e-16 ***
## XWMED  0.021024   0.021150   0.994  0.32054    
## XWFED -0.001782   0.021202  -0.084  0.93305    
## XUN   -0.016026   0.017090  -0.938  0.34870    
## XCIT   0.042237   0.017969   2.351  0.01901 *  
## XAX   -0.044889   0.020399  -2.201  0.02808 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4511 on 735 degrees of freedom
## Multiple R-squared:  0.8011, Adjusted R-squared:  0.7962 
## F-statistic: 164.5 on 18 and 735 DF,  p-value: < 2.2e-16

beta.mle = solve(t(X)%*%X)%*%t(X)%*%y
sig2.mle = mean((y-X%*%beta.mle)^2)
se.mle = sqrt(sig2.mle*diag(solve(t(X)%*%X)))
cbind(beta.mle,se.mle)

##                       se.mle
## LFP  -0.088520910 0.02947983
## WHRS  0.132256030 0.02813058
## KL6   0.019940904 0.01937402
## K618  0.114282051 0.01872898
## WA    0.055324537 0.03766104
## WE   -0.005336829 0.02351607
## WW    0.062273173 0.02477561
## RPWG  0.036699380 0.02505534
## HHRS  0.059382777 0.02077455
## HA    0.041883154 0.03621236
## HE    0.038256653 0.02247579
## HW    0.147562664 0.03157596
## MTR  -0.711872485 0.03249591
## WMED  0.021023693 0.02089564
## WFED -0.001781634 0.02094717
## UN   -0.016025661 0.01688478
## CIT   0.042236942 0.01775312
## AX   -0.044889326 0.02015369

L = beta.mle +qnorm(0.025)*se.mle
U = beta.mle +qnorm(0.975)*se.mle

plot(beta.mle,ylim=range(L,U),xlab="Covariate",ylab="Coefficient",pch=16)
abline(h=0,lty=2)
for (i in 1:p)
  segments(i,L[i],i,U[i])

4 Prior hyperparameters

b0       = rep(0,p)
B0       = diag(1,p)
nu0      = 1
sig20    = 0.2

nu0sig20 = nu0*sig20
iB0      = solve(B0)
iB0b0    = iB0%*%b0

draw = sqrt(1/rgamma(100000,nu0/2,nu0*sig20/2))
quantile(draw,seq(0.1,0.9,by=0.1))

##       10%       20%       30%       40%       50%       60%       70%       80% 
## 0.2716682 0.3485264 0.4309841 0.5330863 0.6656844 0.8586874 1.1702234 1.7785396 
##       90% 
## 3.6269915

5 Posterior sufficient statistics

B1    = solve(iB0+t(X)%*%X)
tcB1  = t(chol(B1))
b1    = B1%*%(iB0b0+t(X)%*%y)
nu1   = nu0 + (n-p)
sig21 = (nu0sig20 + t(y-X%*%b1)%*%y + t(b0-b1)%*%iB0b0)/nu1
se.bayes = sqrt(diag(sig21[1,1]*B1))

L1 = b1+qt(0.025,df=nu1)*se.bayes
U1 = b1+qt(0.975,df=nu1)*se.bayes

par(mfrow=c(1,2))
plot(b1,ylim=range(L,U,L1,U1),xlab="Covariate",ylab="Coefficient",pch=16)
abline(h=0,lty=2)
for (i in 1:p){
  segments(i,L1[i],i,U1[i])
  segments(i+0.2,L[i],i+0.2,U[i],col=2)
  points(i+0.2,beta.mle[i],pch=16,col=2)
}
legend("bottomleft",legend=c("MLE","BAYES"),col=2:1,lwd=2)

sigma2 = seq(0.15,0.25,length=1000)
plot(sigma2,dgamma(1/sigma2,nu1/2,nu1*sig21/2)/(sigma2^2),xlab=expression(sigma^2),ylab="Density",type="l")
lines(sigma2,dgamma(1/sigma2,nu0/2,nu0*sig20/2)/(sigma2^2),lty=2)
abline(v=sig2.mle,col=2,lwd=2)
legend("topleft",legend=c("Prior","Posterior","MLE"),col=c(1,1,2),lty=c(1,2,1),lwd=2)

6 Predictive - full model (M0)

m   = X%*%b0
V   = sig20*(diag(1,n)+X%*%B0%*%t(X))
res = t(chol(V))%*%(y-m)
predictive0 = sum(dt(res,df=nu0,log=TRUE))
predictive0

## [1] -2253.028

7 Predictive - reduced model (M1)

X1 = X[,c(1,2,4,7,9,12,13,17,18)]
p1 = ncol(X1)
m   = X1%*%b0[1:p1]
V   = sig20*(diag(1,n)+X1%*%B0[1:p1,1:p1]%*%t(X1))
res = t(chol(V))%*%(y-m)
predictive1 = sum(dt(res,df=nu0,log=TRUE))
predictive1

## [1] -2107.195

8 Predictive - another reduced model (M2)

X2 = X[,c(1,2,4,12,13)]
p2 = ncol(X2)
m   = X2%*%b0[1:p2]
V   = sig20*(diag(1,n)+X2%*%B0[1:p2,1:p2]%*%t(X2))
res = t(chol(V))%*%(y-m)
predictive2 = sum(dt(res,df=nu0,log=TRUE))
predictive2

## [1] -2015.077

9 Computing log Bayes factors

logB10 = predictive1-predictive0
logB20 = predictive2-predictive0
logB21 = predictive2-predictive1

c(predictive0,predictive1,predictive2)

## [1] -2253.028 -2107.195 -2015.077

c(logB10,logB20,logB21)

## [1] 145.83288 237.95107  92.11819

exp(c(logB10,logB20,logB21))

## [1]  2.159812e+63 2.191989e+103  1.014898e+40

10 Final model

X2       = X[,c(1,2,4,12,13)]
names1   = names[c(1,2,4,12,13)]
p        = 5
b0       = rep(0,p)
B0       = diag(1,p)
nu0      = 1
sig20    = 0.2
nu0sig20 = nu0*sig20
iB0      = solve(B0)
iB0b0    = iB0%*%b0

B1    = solve(iB0+t(X2)%*%X2)
tcB1  = t(chol(B1))
b1    = B1%*%(iB0b0+t(X2)%*%y)
nu1   = nu0 + (n-p)
sig21 = (nu0sig20 + t(y-X2%*%b1)%*%y + t(b0-b1)%*%iB0b0)/nu1
se.bayes = sqrt(diag(sig21[1,1]*B1))

L1 = b1+qt(0.025,df=nu1)*se.bayes
U1 = b1+qt(0.975,df=nu1)*se.bayes

par(mfrow=c(1,1))
plot(b1,ylim=range(L1,U1),xlab="Covariate",ylab="Coefficient",pch=16,axes=FALSE)
axis(2);box()
axis(1,at=1:p,lab=names1)
abline(h=0,lty=2)
for (i in 1:p)
  segments(i,L1[i],i,U1[i])
title("Final Bayesian linear model")

Regressors:

• LFP “A dummy variable = 1 if woman worked in 1975, else 0”;

• WHRS “Wife’s hours of work in 1975”;

• K618 “Number of children between ages 6 and 18 in household”;

• HW “Husband’s wage, in 1975 dollars”;

• MTR “This is the marginal tax rate facing the wife, and is taken from published federal tax tables (state and local income taxes are excluded). The taxable income on which this tax rate is calculated includes Social Security, if applicable to wife.”;

11 Residual analysis

fitted   = X2%*%b1
residual = (y-fitted)/sqrt(sig21[1,1])
par(mfrow=c(1,2))
plot(residual,xlab="Observation",ylab="Standardized residuals")
abline(h=2,lty=2)
abline(h=0,lty=2)
abline(h=-2,lty=2)
title(paste("Residuals outside [-2,2] = ",round(100*mean(abs(residual)>2),1),"%",sep=""))

plot(y,fitted,xlim=range(y,fitted),ylim=range(y,fitted),xlab="Observed response",ylab="Fitted response")
abline(0,1,col=2,lwd=2)