Simple linear regression

1 Weight vs height

We collected data from 15 volunteers (my students and myself) and below we summarize the data information graphically and entertain two competing simple linear regression models: i) Gaussian errors, ii) Student’s \(t\) errors.

hedi = 115
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,hedi),n,2,byrow=TRUE)

height = data[,1]
weight = data[,2]

par(mfrow=c(1,1))
plot(height,weight,xlab="Height (in cm)",ylab="Weight (in kg)")
points(height[n],weight[n],col=2,pch=16)
abline(lm(weight[1:(n-1)]~height[1:(n-1)]),lwd=2)
abline(lm(weight~height),col=2,lwd=2)

1.1 Fitted lines

Below we use the lm (linear model) function from R to run MLE/OLS for the regression \[ \mbox{weight}_i = \alpha + \beta \ \mbox{height}_i + \varepsilon_i \qquad\qquad \varepsilon_i \ \ \mbox{iid} \ \ N(0,\sigma^2). \] The two fitted lines, including and excluing the outlier (observation 15) are presented in the above picture.

Let us explore how the outlier affects the linear regression fits.

• The adjusted \(R^2\) goes from 33% to 43% (including outlier).

• The height coefficient is estimated at \(0.52\) with p-value at \(0.02\) for the null hypothesis \(H_0:\beta=0\), when excluding the outlier.

• The height coefficient is estimated at \(0.84\) with p-value at \(0.005\) for the null hypothesis \(H_0:\beta=0\), when including the outlier.

• Therefore, at the 1% level, one would rejected the null hypothesis when the outlier is included and fail to reject when the oulier is excluded.

• The estimates of \(\sigma\) are \(6.7\) (excluding the outlier) and \(9.7\) (including the outlier), ie. a significant increase in the variability of the unexplained error when the outlier is included.

summary(lm(weight[1:(n-1)]~height[1:(n-1)]))

## 
## Call:
## lm(formula = weight[1:(n - 1)] ~ height[1:(n - 1)])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.889  -5.197   1.994   4.308  10.011 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       -16.0385    34.2319  -0.469   0.6478  
## height[1:(n - 1)]   0.5168     0.1915   2.699   0.0194 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.726 on 12 degrees of freedom
## Multiple R-squared:  0.3777, Adjusted R-squared:  0.3258 
## F-statistic: 7.283 on 1 and 12 DF,  p-value: 0.01936

summary(lm(weight~height))

## 
## Call:
## lm(formula = weight ~ height)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.262  -7.108  -1.589   5.960  22.536 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -71.5004    44.5970  -1.603  0.13289   
## height        0.8366     0.2478   3.375  0.00497 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.66 on 13 degrees of freedom
## Multiple R-squared:  0.4671, Adjusted R-squared:  0.4261 
## F-statistic: 11.39 on 1 and 13 DF,  p-value: 0.004971

2 Bayesian inference

loglike.g = function(alpha,beta,sig,x,y){
  sum(dnorm(y,alpha+beta*x,sig,log=TRUE))
}
loglike.t = function(alpha,beta,sig,x,y){
    tau = sqrt((nu-2)*sig^2/nu)
  sum(dt((y-alpha-beta*x)/tau,df=nu,log=TRUE))-n*log(tau)
}


y  = weight
x  = height
y1 = y[1:(n-1)]
x1 = x[1:(n-1)]

nu      = 3
M       = 100000
alpha.d = runif(M,-250,150)
beta.d  = runif(M,0,2)
sig.d   = runif(M,1,20)
w       = matrix(0,M,2)
w1      = matrix(0,M,2)
for (i in 1:M){
  w[i,1] = loglike.g(alpha.d[i],beta.d[i],sig.d[i],x,y)
  w[i,2] = loglike.t(alpha.d[i],beta.d[i],sig.d[i],x,y)
  w1[i,1] = loglike.g(alpha.d[i],beta.d[i],sig.d[i],x1,y1)
  w1[i,2] = loglike.t(alpha.d[i],beta.d[i],sig.d[i],x1,y1)
}
ww  = exp(w)
ww1 = exp(w1)
w[,1]   = exp(w[,1]-max(w[,1]))
w[,2]   = exp(w[,2]-max(w[,2]))
w1[,1]  = exp(w1[,1]-max(w1[,1]))
w1[,2]  = exp(w1[,2]-max(w1[,2]))
ind.n   = sample(1:M,size=M,replace=TRUE,prob=w[,1])
ind.t   = sample(1:M,size=M,replace=TRUE,prob=w[,2])
ind1.n  = sample(1:M,size=M,replace=TRUE,prob=w1[,1])
ind1.t  = sample(1:M,size=M,replace=TRUE,prob=w1[,2])

pars = array(0,c(2,M,3))
pars[1,,] = cbind(alpha.d[ind.n],beta.d[ind.n],sig.d[ind.n])
pars[2,,] = cbind(alpha.d[ind.t],beta.d[ind.t],sig.d[ind.t])

pars1 = array(0,c(2,M,3))
pars1[1,,] = cbind(alpha.d[ind1.n],beta.d[ind1.n],sig.d[ind1.n])
pars1[2,,] = cbind(alpha.d[ind1.t],beta.d[ind1.t],sig.d[ind1.t])

names = c("alpha","beta","sigma")
model = c("Gaussian errors","Student's t errors")

2.1 Marginal posterior distributions (with outlier)

par(mfrow=c(2,3))
for (j in 1:2)
  for (i in 1:3)
    hist(pars[j,,i],prob=TRUE,main=names[i],xlab=model[j])

2.2 Marginal posterior distributions (without outlier)

par(mfrow=c(2,3))
for (j in 1:2)
  for (i in 1:3)
    hist(pars1[j,,i],prob=TRUE,main=names[i],xlab=model[j])

2.3 Comparison

par(mfrow=c(2,3))
for (j in 1:2)
  for (i in 1:3)
    boxplot(pars[j,,i],pars1[j,,i],main=names[i],xlab=model[j],outline=FALSE,
            ylim=range(pars[,,i],pars1[,,i]),names=c("with outlier","without outlier"))

par(mfrow=c(3,1))
for (i in 1:3)
    boxplot(pars[1,,i],pars1[1,,i],pars[2,,i],pars1[2,,i],main=names[i],outline=FALSE,
            ylim=range(pars[,,i],pars1[,,i]),
            names=c("G with outlier","G without outlier","t with outlier","t without outlier"))

2.4 Predictive, Bayes factor and posterior model probability

2.4.1 With outlier

The Bayes factor is around 0.57, which indicates that both Gaussian and Student’t quite similar in terms of model adequacy. The posterior model probabilities are roughly 1/3 for the Gaussian model and 2/3 for the Student’s \(t\) model.

pred  = apply(ww,2,mean)
pmp   = pred/sum(pred)
B12   = pred[1]/pred[2]

pred

## [1] 3.191276e-27 5.511371e-27

B12

## [1] 0.5790348

pmp

## [1] 0.3667018 0.6332982

2.4.2 Without outlier

When removing the outlier, the Bayes factor is around 10, which indicates that the Gaussian model is moderately better than theStudent’t model. The posterior model probabilities are roughly 90% for the Gaussian model and 10% for the Student’s \(t\) model. In words, the presence of the outlier gives the Student’s \(t\) model a fitting strengh. Such strengh disappear when the outlier is removed and the Gaussian model suffices as a linear model.

pred1 = apply(ww1,2,mean)
pmp1  = pred1/sum(pred1)
B121  = pred1[1]/pred1[2]

pred1

## [1] 7.521713e-24 7.705170e-25

B121

## [1] 9.761904

pmp1

## [1] 0.90707964 0.09292036