Simple linear regression with Gaussian errors

1 Weight and height dataset

We collected data from 15 volunteers (my students and myself).

n    = 15
data = matrix(c(178,70,175,80,184,82,200,89,183,80,173,65,185,70,162,70,
190,85,174,80,180,87,165,74,176,72,174,63,196,113),n,2,byrow=TRUE)
height = scale(data[,1])
weight = scale(data[,2])
par(mfrow=c(1,1))
plot(height,weight,xlab="Height (standardized)",ylab="Weight (standardized)")
points(height[n],weight[n],col=2,pch=16,cex=2)
abline(v=0,lty=2)
abline(h=0,lty=2)

2 Maximum likelihood estimation (MLE)

Below we use the lm (linear model) function from R to run the OLS estimation for the regression \[ y_i = \beta \ x_i + \varepsilon_i \qquad\qquad \varepsilon_i \ \ \mbox{iid} \ \ N(0,\sigma^2), \] where \(x_i=\mbox{height}_i\) and \(y_i=\mbox{weight}_i\), for \(i=1,\ldots,n\). From the frequentist/Fisherian approach, we have learned that the MLEs of \(\beta\) and \(\sigma^2\) are, respectively, \[ {\widehat \beta} = \frac{\sum_{i=1}^n y_i x_i}{\sum_{i=1}^n x^2_i} \qquad \mbox{and} \qquad {\widehat \sigma}^2 = \frac{\sum_{i=1}^n (y_i-{\widehat y}_i)^2}{n}. \] We also know the sampling distribution of these estimators: \[ {\widehat \beta}|{\widehat \sigma}^2 \sim t_{n-1}\left(\beta,\frac{{\widehat \sigma}^2}{\sum_{i=1}^n x_i^2}\right) \qquad and \qquad {\widehat \sigma}^2|\sigma^2 \sim IG\left(\frac{n-1}{2},\frac{(n-1)\sigma^2}{2}\right), \] from which confidence intervals and hypothesis testing can be derived. An unbiased estimator of \(\sigma^2\) is \(s^2 = n{\widehat \sigma}^2/(n-1)\).

bhat  = sum(height*weight)/sum(height^2)
sighat = sqrt(sum((weight-bhat*height)^2)/n)
s = sqrt(n*sighat^2/(n-1))
se = s/sqrt(sum(height^2))
round(c(bhat,se),4)

## [1] 0.6875 0.1941

round(c(sighat,s),4)

## [1] 0.7016 0.7262

# The -1 is to inform the lm function that the regression has no intercept
summary(lm(weight~height-1))

## 
## Call:
## lm(formula = weight ~ height - 1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0537 -0.5485 -0.1120  0.4731  1.7021 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)   
## height   0.6875     0.1941   3.542  0.00325 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7262 on 14 degrees of freedom
## Multiple R-squared:  0.4727, Adjusted R-squared:  0.435 
## F-statistic: 12.55 on 1 and 14 DF,  p-value: 0.003249

2.1 OLS fit with and without the (red) outlier

par(mfrow=c(1,1))
plot(height,weight,xlab="Height (standardized)",ylab="Weight (standardized)")
points(height[n],weight[n],col=2,pch=16,cex=2)
abline(lm(weight~height-1),col=2,lwd=2)
abline(lm(weight[1:(n-1)]~height[1:(n-1)]-1),lwd=2)
legend("topleft",legend=c("With outlier","Without outlier"),col=2:1,lwd=2,lty=1,bty="n")
abline(v=0,lty=2)
abline(h=0,lty=2)

3 Bayesian inference

Let us assume the following independent prior for \(\beta\) and \(\sigma^2\) \[ p(\beta,\sigma^2) = p(\beta|\sigma^2)p(\sigma^2), \] with \(\beta \sim N(b_0,\sigma^2B_0)\) and \(\sigma^2 \sim IG(c_0,d_0)\), with hyperparameters \(b_0=0\), \(B_0=1\), \(c_0=2\) and \(d0=1\). Notice that this is a conjugate prior for \((\beta,\sigma^2)\), so closed-form posterior inference is readily available.

3.1 Draws from the prior

b0=0
B0=1
c0=2
d0=1
M    = 1000
sig2 = 1/rgamma(M,c0,d0)
beta = rnorm(M,b0,sqrt(sig2*B0))
par(mfrow=c(1,1))
boxplot(cbind(beta,sqrt(sig2)),names=c(expression(beta),expression(sigma)),outline=F,xlab="",horizontal=TRUE)
points(bhat,1,col=2,pch=16,cex=2)
points(s,2,col=2,pch=16,cex=2)
title("Prior draws (MLE in red)")

4 Your homework is listed below

0. Derive the closed form posterior distribution \(p(\beta,\sigma^2|\mbox{data})\), where \(\mbox{data}=\{(x_1,y_1),\ldots,(x_n,y_n)\}\).

1. Implement the sampling importance resampling (SIR) algorithm to obtain draws from the posterior \(p(\beta,\sigma^2|\mbox{data})\). You can either use the prior for \((beta,\sigma^2)\) as proposal density or you can use the following proposal: \(q(\beta) \equiv N(0.7,0.25)\) and \(q(\sigma^2) \equiv IG(7,5)\).

• Note: These proposal densities are inspired by the sampling distributions of the corresponding estimators. Recall that the proposal distributions can be anything with or without the use of the data since it is not altering the prior in any way.

2. Implement the Gibbs sampler to obtain draws from the posterior \(p(\beta,\sigma^2|\mbox{data})\). Check the mixing of the Markov chains and the convergence of the draws to the posterior.

3. Repeat 2) but replace the Gibbs sampler by a Metropolis algorithm with proposals \(q(\beta) \equiv N(0.7,0.25)\) and \(q(\sigma^2) \equiv IG(7,5)\), which already appear in 1).

4. If the three algorithms are running properly, the results should be the same up to small Monte Carlo errors.