Exercise

  1. Use $\beta_0=1$, $\sigma^2_{\beta}=0.1$ and $\sigma^2_e = 1.0$ to generate a data set with 10 observations from model (2) with $k=15$ covariates.

  2. Setup and solve the mixed model equations given by (3).

  3. Sample the elements of $\boldsymbol{\beta}$ using Gibbs.

  4. Compute the posterior mean of $\boldsymbol{\beta}$ from the samples and compare with the mixed model solutions.

  5. Compute the posterior covariance matrix from the sampled values. Compare results with inverse of the mixed-model coefficient matrix.

Simulation of Data

In [12]:
using Distributions
using StatsBase
In [13]:
n = 10 #number of observations
k = 15    #number of covariates

varBeta = 0.1
varRes  = 1.0

x = sample([0,1,2],(n,k))
X = [ones(n) x]

betaTrue=[1, randn(k)*sqrt(varBeta)]
y = X*betaTrue+ randn(n)*sqrt(varRes);

MME

In [14]:
λ=varRes/varBeta
d=eye(k+1)
d[1,1]=0
lhs=X'X+d*λ
rhs=X'y
ilhs = inv(lhs)
sol=ilhs*rhs
Out[14]:
16-element Array{Float64,1}:
  0.739611 
 -0.0985318
 -0.0733791
  0.244959 
  0.171841 
  0.0126194
  0.0451185
 -0.0675262
 -0.0307764
 -0.181414 
 -0.193714 
  0.255449 
 -0.194699 
 -0.183016 
 -0.113338 
  0.153745 
In [15]:
ilhs*varRes
Out[15]:
16x16 Array{Float64,2}:
  1.477      -0.0714228    -0.117995     …  -0.101212     -0.05469   
 -0.0714228   0.0758958     0.00701886       0.00749083    0.0067372 
 -0.117995    0.00701886    0.0702412        0.000819962  -0.00193193
 -0.0733976   0.00211372    0.00179398       0.00259769   -0.00573441
 -0.0891733   0.00406518    0.00509045       0.000641145  -0.00204173
 -0.127917    0.000686534   0.0153226    …   0.00329912   -0.0067671 
 -0.0275212  -0.000976669   0.0147982       -0.0155726     0.00236965
 -0.107998   -0.00026764    0.0087143        0.006673      0.00504872
 -0.0340253  -0.00482864    0.00436116       0.00322588    0.00317027
 -0.0642605  -0.0210254     0.0034216       -0.000252368  -0.00144475
 -0.0954171   0.00105101    0.00223758   …  -0.00765011    0.0131453 
 -0.126453    0.00377223   -0.00606176       0.0139139    -0.00279948
 -0.0945257  -0.00416139   -0.00977513      -0.00482177    0.00724818
 -0.0921117   0.000430086  -0.00751776      -0.00323865    0.00813793
 -0.101212    0.00749083    0.000819962      0.073409     -0.00654971
 -0.05469     0.0067372    -0.00193193   …  -0.00654971    0.0779866 

MCMC:

In [19]:
# loop for Gibbs sampler
niter  = 100000 # number of samples
burnin = 1000
b     = zeros(1+k)
meanB = zeros(1+k)
vB    = zeros(1+k,1+k)
meansigma2 = 0
nu = 4
sigma2 = 2
covariates = length(b) - 1 #Number of covariates
b_mcmc = zeros(niter, length(b))

for iter = 1:niter
    
    # sampling intercept
    w = y - X*b + X[:,1]*b[1]
    x = X[:,1]
    xpxi = 1/(x'x)[1]
    bHat = (xpxi*x'w)[1] 
    b[1] = rand(Normal(bHat, sqrt(xpxi))) # using residual var = 1 

    #sample covariates
    for j = 2:length(b)
        w = y - X*b + X[:,j]*b[j]
        x = X[:,j]
        xpxi = 1/((x'x)[1]+λ)
        bHat = (xpxi*x'w)[1] 
        b[j] = rand(Normal(bHat, sqrt(xpxi)))
    end
    if iter >= burnin
        b_mcmc[iter,:] = b
        meanB += b
        vB    += b*b'
    end
end

Posterior Mean of $\beta$ Estimated from MCMC Samples

In [20]:
meanB /= (niter-burnin)
Out[20]:
16-element Array{Float64,1}:
  1.20281   
 -0.0941883 
 -0.143509  
  0.314295  
  0.222457  
 -0.108558  
  0.0911845 
 -0.149743  
 -0.00393558
 -0.263458  
 -0.296184  
  0.348927  
 -0.210706  
 -0.341597  
 -0.152182  
  0.179394  

Posterior Covariance of $\beta$ Estimated from MCMC Samples

In [9]:
vB/(niter-burnin) - meanB*meanB'
Out[9]:
16x16 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0

Unknown variance components

In [29]:
k=15
niter  = 100000 # number of samples
burnin = 1000
b     = zeros(1+k)
meanB = zeros(1+k)
vB    = zeros(1+k,1+k)
meansigma2 = 0
nu = 4
sigma2 = 2
covariates = length(b) - 1 #Number of covariates

dfEffectVar=4
scaleVar=0.1
chi2=Chisq(dfEffectVar+covariates)

sampleVar = zeros(niter)

for iter = 1:niter
    
    # sampling intercept
    w = y - X*b + X[:,1]*b[1]
    x = X[:,1]
    xpxi = 1/(x'x)[1]
    bHat = (xpxi*x'w)[1] 
    b[1] = rand(Normal(bHat, sqrt(xpxi))) # using residual var = 1 

    #sample covariates
    λ=varRes/varBeta
    for j = 2:length(b)
        w = y - X*b + X[:,j]*b[j]
        x = X[:,j]
        xpxi = 1/((x'x)[1]+λ)
        bHat = (xpxi*x'w)[1] 
        b[j] = rand(Normal(bHat, sqrt(xpxi)))
    end
    
    varBeta = (scaleVar*dfEffectVar + dot(b[2:end],b[2:end]))/rand(chi2)  

    if iter >= burnin
        meanB += b
        vB    += b*b'
    end
    
    sampleVar[iter]=varBeta
end
In [34]:
using Gadfly
plot(x=sampleVar,Geom.histogram,Guide.title("Posterior distribution of varBeta"),
Guide.ylabel("Frequency"),
Guide.xlabel("varBeta"))
Out[34]:
varBeta -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -1.50 -1.45 -1.40 -1.35 -1.30 -1.25 -1.20 -1.15 -1.10 -1.05 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 2.35 2.40 2.45 2.50 2.55 2.60 2.65 2.70 2.75 2.80 2.85 2.90 2.95 3.00 -2 0 2 4 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -1.5×10⁴ -1.0×10⁴ -5.0×10³ 0 5.0×10³ 1.0×10⁴ 1.5×10⁴ 2.0×10⁴ 2.5×10⁴ -1.00×10⁴ -9.50×10³ -9.00×10³ -8.50×10³ -8.00×10³ -7.50×10³ -7.00×10³ -6.50×10³ -6.00×10³ -5.50×10³ -5.00×10³ -4.50×10³ -4.00×10³ -3.50×10³ -3.00×10³ -2.50×10³ -2.00×10³ -1.50×10³ -1.00×10³ -5.00×10² 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ 1.25×10⁴ 1.30×10⁴ 1.35×10⁴ 1.40×10⁴ 1.45×10⁴ 1.50×10⁴ 1.55×10⁴ 1.60×10⁴ 1.65×10⁴ 1.70×10⁴ 1.75×10⁴ 1.80×10⁴ 1.85×10⁴ 1.90×10⁴ 1.95×10⁴ 2.00×10⁴ -1×10⁴ 0 1×10⁴ 2×10⁴ -1.0×10⁴ -9.0×10³ -8.0×10³ -7.0×10³ -6.0×10³ -5.0×10³ -4.0×10³ -3.0×10³ -2.0×10³ -1.0×10³ 0 1.0×10³ 2.0×10³ 3.0×10³ 4.0×10³ 5.0×10³ 6.0×10³ 7.0×10³ 8.0×10³ 9.0×10³ 1.0×10⁴ 1.1×10⁴ 1.2×10⁴ 1.3×10⁴ 1.4×10⁴ 1.5×10⁴ 1.6×10⁴ 1.7×10⁴ 1.8×10⁴ 1.9×10⁴ 2.0×10⁴ Frequency Posterior distribution of varBeta