Applied Statistical Methods - Notebook 6

Problem 1: Example with LASSO

The file available at https://charlotte-ngs.github.io/gelasmss2021/data/asm_ex06_p01_lasso.txt contains a dataset with genotypes from 100 SNP-Loci. In addition to the genomic information, the dataset also holds observations of a certain trait for 50 animals. The dataset can be read into a matrix in R with the following statement.

mat_lasso_data <- matrix(scan("https://charlotte-ngs.github.io/gelasmss2021/data/asm_ex06_p01_lasso.txt"), nrow = 50, byrow = TRUE)

Let us have a look at the first 5 rows and the first 5 columns of the matrix that stores the dataset.

mat_lasso_data[1:n_nr_row,1:n_nr_col]

          [,1] [,2] [,3] [,4] [,5]
[1,] -40.39872   -1    0   -1   -1
[2,] -46.35871   -1   -1   -1    0
[3,] -33.60278   -1   -1   -1   -1
[4,] -48.47177    0   -1   -1   -1
[5,] -38.82089   -1   -1    0   -1

From this output, we can see that the observations of all animals can be found in the first column of the data matrix mat_lasso_data. In columns 2 to 101 of the data matrix there are the genotypes of all SNP loci. The linear model is fitted with LASSO using the function glmnet() from the package glmnet. The SNP genotypes are used as explanatory variables and the observations are the response variables.

Your Tasks

• Use the following R-Statement to estimate the SNP-effects using LASSO

require(glmnet)
fitsnp <- glmnet(x = mat_lasso_data[, -1], y = mat_lasso_data[, 1])

• Visualize the dependency between the value of the penalty term \(\lambda\) and the number of explanatory variables which are not \(0\).

plot(fitsnp, xvar = "lambda", label = TRUE)

• Use a cross-validation to determine the value of \(\lambda\).

cvfitsnp <- cv.glmnet(x = mat_lasso_data[, -1], y = mat_lasso_data[, 1])

• Show the results of the cross-validation in a plot using the function plot().

plot(cvfitsnp)

• In the plot of the cross-validation results there are two dashed lines which both indicated special values for \(\lambda\). The first value is the minimum of all \(\lambda\)-values and the second is the one that sets the most explanatory variables to \(0\) with the restriction that the sum of squared errors is not further away than one standard deviation from its minimum. The two \(\lambda\)-values are obtained with

cvfitsnp$lambda.min
cvfitsnp$lambda.1se

• Find all coefficients which are not \(0\) for both \(\lambda\)-values and compare them to the true values taken from the simulation.

coefmin <- coef(cvfitsnp, s = "lambda.min")
(cofminnz <- coefmin[coefmin[, 1] != 0,])

coef1se <- coef(cvfitsnp, s = "lambda.1se")
(coef1senz <- coef1se[coef1se[, 1] != 0, ])

The true SNP-positions from the simulation are:

(vec_sign_snp_idx <- c(73,54,26,30,7))

[1] 73 54 26 30  7

Your Solution

Run the R-statements as they are shown in the hints. This will result in a set of SNP-positions with an effect that is different from \(0\). These SNP-positions can then be compared to the true positions from the simulation. The LASSO-estimates of the SNP-effects are obtained with the following steps.

# fit the linear model with glmnet

# Visualize the results using the plot() function

# Determine lambda using cross-validation with the function cv.glmnet()

# Show the results of the cross-validation using the plot() function

# determine the two lambda values, given in the above plot

# determine the non-zero coefficients for the SNP-effects which are not zero for the two lambda values

# compare the SNP-positions with a non-zero effect to the true SNP positions used in the simulation

Problem 2: Bayesian Regression Analysis

Given is the earlier used dataset of breast circumference and body weight.

Dataset for Regression of Body Weight on Breast Circumference for ten Animals

Animal	Breast Circumference	Body Weight
1	176	471
2	177	463
3	178	481
4	179	470
5	179	496
6	180	491
7	181	518
8	182	511
9	183	510
10	184	541

The model that is used is a simple linear regression model given by

\[y_i = \beta_0 + \beta_1 * x_i + \epsilon_i\].

where \(y_i\) corresponds to the body weight of animal \(i\), \(x_i\) is the breast circumference of animal \(i\), \(\beta_0\) is the unknown intercept and \(\beta_1\) is the unknown regression coefficient. For reasons of simplicity, we assume the residual variance \(\sigma^2\) to be known. For the later computations, we insert the estimate that is obtained from the lm() function. This value corresponds to \(\sigma^2 = 122.8\).

Bayesian Estimation Of Unknowns

As already mentioned during the lecture, Bayesian estimates of unknowns are based on the posterior distribution of the unknowns given the knowns. For our regression model the unknowns correspond to

\[\beta = \left[\begin{array}{c}\beta_0 \\ \beta_1 \end{array}\right]\]

The posterior distribution of the unknowns given the knowns is \(f(\beta | y)\). Using Bayes’ Theorem we can write \(f(\beta | y)\) as

\[\begin{align} f(\beta | y) & = \frac{f(\beta, y)}{f(y)} \notag \\ & = \frac{f(y | \beta)f(\beta)}{f(y)} \notag \\ & \propto f(y | \beta)f(\beta) \notag \end{align}\]

When we do not have any specific prior knowledge about \(\beta\), the prior distribution \(f(\beta)\) for the unknown \(\beta\) is set to a constant. Therefore we can write

\[\begin{align} f(\beta | y) & \propto f(y | \beta)f(\beta) \notag \\ & \propto f(y | \beta) \notag \end{align}\]

Assuming a normal distribution for the data causes the likelihood \(f(y | \beta)\) to be a multivariate normal distribution.

\[\begin{align} f(\beta | y) & \propto f(y | \beta) \notag \\ & = (2\pi\sigma^2)^{-n/2} exp \left\{ -{1 \over 2} \frac{(y - X\beta)^T(y - X\beta)}{\sigma^2} \right\} (\#eq:BayesLikelihood) \end{align}\]

The above expression @ref(eq:BayesLikelihood) is an \(n-\) dimensional normal distribution with expected value \(X\beta\) and variance-covariance matrix corresponding to \(I\sigma^2\). But because we have just two unknowns \(\beta_0\) and \(\beta_1\) the posterior distribution \(f(\beta | y)\) must have two dimensions and not \(n\). The following re-arrangement can solve this problem. Let us set the variable \(Q\) to

\[Q = (y - X\beta)^T(y - X\beta) = y^Ty - 2y^TX\beta + \beta^T(X^TX)\beta\]

Introducing the least squares estimate \(\hat{\beta} = (X^TX)^{-1}X^Ty\) into the above equation by replacing \(y^TX\) with \(\hat{\beta}^T(X^TX)\) results in

\[Q = y^Ty - 2 \hat{\beta}^T(X^TX)\beta + \beta^T(X^TX)\beta = y^Ty + (\beta-\hat{\beta})^T(X^TX)(\beta-\hat{\beta}) - \hat{\beta}^T(X^TX)\hat{\beta}\]

Inserting this last result back into @ref(eq:BayesLikelihood) gives

\[\begin{align} f(\beta | y) & \propto f(y | \beta) \notag \\ & = (2\pi\sigma^2)^{-n/2} exp \left\{ -{1 \over 2} \frac{(y - X\beta)^T(y - X\beta)}{\sigma^2} \right\} \notag \\ & = (2\pi\sigma^2)^{-n/2} exp \left\{ -{1 \over 2} \frac{y^Ty + (\beta-\hat{\beta})^T(X^TX)(\beta-\hat{\beta}) - \hat{\beta}^T(X^TX)\hat{\beta}}{\sigma^2} \right\} \notag \\ & = (2\pi\sigma^2)^{-n/2} \left[exp \left\{ -{1 \over 2} \frac{y^Ty}{\sigma^2}\right\} * exp\left\{ -{1 \over 2} \frac{ (\beta-\hat{\beta})^T(X^TX)(\beta-\hat{\beta})}{\sigma^2} \right\} * exp\left\{ -{1 \over 2} \frac{ - \hat{\beta}^T(X^TX)\hat{\beta}}{\sigma^2} \right\} \right] \notag \\ & \propto exp\left\{ -{1 \over 2} \frac{ (\beta-\hat{\beta})^T(X^TX)(\beta-\hat{\beta})}{\sigma^2} \right\} (\#eq:BayesLikelihoodRef) \end{align}\]

The last proportionality results from the fact that only the term depending on \(\beta\) is retained. All other terms not depending on \(\beta\) are constant factors with respect to \(\beta\) and can therefore be dropped. Thus \(f(\beta|y)\) can be written as

\[f(\beta | y) \propto exp\left\{ -{1 \over 2} \frac{ (\beta-\hat{\beta})^T(X^TX)(\beta-\hat{\beta})}{\sigma^2} \right\}\] which is recognized as proportional to a two dimensional normal density with mean \(\hat{\beta}\) and variance \((X^TX)^{-1}\sigma^2\). Thus in the simple setting the mean of the posterior mean can already be seen from the above formula. But in a more complex setting, the posterior distribution does not have a standard form and we need to setup a sampling scheme which allows us to draw random numbers from the posterior distribution. The sampling scheme that we are introducing here is called the Gibbs Sampler.

Gibbs Sampler for \(\beta\)

The simple regression model that we are using for the breast circumference and the body weight data can be written in matrix-vector notation as

\[y = 1\beta_0 + x\beta_1 + \epsilon\]

In the Gibbs sampling scheme both unknowns \(\beta_0\) and \(\beta_1\) are sampled from their full conditional distributions. For \(\beta_0\) the full conditional posterior distribution is \(f(\beta_0 | \beta_1, y)\) which is computed for the current value of \(\beta_1\). Separating \(\beta_0\) from the other unknowns yields the linear model

\[w_0 = 1\beta_0 + \epsilon\] where \(w_0 = y - x\beta_1\). The least squares estimator of \(\beta_0\) is

\[\hat{\beta}_0 = (1^T1)^{-1}1^Tw_0\] with variance

\[var(\hat{\beta}_0) = (1^T1)^{-1} \sigma^2\]

Applying the same strategy as for \(f(\beta | y)\), it can be shown that \(f(\beta_0 | \beta_1, y)\) is a normal distribution with mean \(\hat{\beta}_0\) as mean and \((1^T1)^{-1} \sigma^2\) as variance. The full-conditional posterior of \(\beta_1\) can be derived the same way, leading to

\[\hat{\beta}_1 = (x^Tx)^{-1}x^Tw_1\]

with variance \(var(\hat{\beta}_1) = (x^Tx)^{-1} \sigma^2\) where \(w_1 = y - 1\beta_0\).

Your Task

• Create a Gibbs Sampling scheme for the dataset shown in Table @ref(tab:dataregression).
• Use the mean of the generated samples as an estimate for the unknowns \(\beta_0\) and \(\beta_1\).

Your Solution

In the problem description, the residual variance \(\sigma^2\) is assumed to be known. Instead of using a fixed value, we use the estimate that can be computed based on the residuals from a least squares model.

# compute residual variance based on residuals of fitting lm() to data

# prepare the vector y of observations and the matrix X for the intercept and the predictor variable

# initialise a vector beta for the current estimate and a vector meanBeta for the mean of all the samples

# loop over a given number of iterations to produce the samples for the unknowns

---

Latest Changes: 2021-04-12 01:02:11 (peter)