Simple Linear Regression

• Ref: https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/06/lecture-06.pdf
• Conditional pdf of \(Y\) for each \(x\) is \(p(y|X = x; \beta_0, \beta_1, \sigma^2)\)

• Probability density under the model \[\prod_{i=1}^{n} p(y_i|x_i;\beta_0,\beta_1,\sigma^2) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi \sigma^2}} e ^{-\frac{(y_i - (\beta_0 + \beta_1 x_i))^2}{2\sigma^2}}\]

• Likelihood given parameters \(b_0, b_1, s^2\) \[\prod_{i=1}^{n} p(y_i|x_i;b_0,b_1,s^2) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi s^2}} e ^{-\frac{(y_i - (b_0 + b_1 x_i))^2}{2s^2}}\]

• Log-likelihood \[L(b_0,b_1,s^2) = -\frac{n}{2}log2\pi - nlogs - \frac{1}{2s^2}\sum_{n=1}^n(y_i - (b_0+b_1x_i))^2\]

Simulate Data

n <- 30
set.seed(123456)
x1 <- rnorm(n = n, mean = 1, sd = 2)
x2 <- rnorm(n = n, mean = 3, sd = 4)
e <- rnorm(n = n, mean = 0, sd = 6)
w <- runif(n = n, min = 0, max = 1)
y <- 3 + x1 * 2 + x2 * 0.5 + e


dt <- data.table(y, x1, x2, w)

X <- as.matrix(cbind(x0 = rep(1, times = n), x1, x2))

p <- dim(X)[2]


dim(X) %>% as.data.table %>% prt(caption = "Dimensions of Design Matrix (X)")

Y <- as.matrix(dt$y)

m <- lm(y ~ x1 + x2, data = dt)
ms <- summary(m)

Matrix Algebra

• Sum of squares

t(X) %*% X %>% prt(caption = "Covariance of Matrix X")

Estimate Beta

\[y = X \beta\] \[X^T y = (X^T X) \beta\] \[\beta = (X^TX)^{-1}X^Ty\]

beta <- solve(t(X) %*% X) %*% t(X) %*% y

cbind(beta, coef(m)) %>%
    prt(caption = "Betas", col.names = c("Matrix Algebra", "lm"))

Estimate Variance (Sigma)

• Ref: https://www.stat.purdue.edu/~boli/stat512/lectures/topic3.pdf
• \(Y\) is a multivariate normal distribution

• Estimate \(\sigma^2\) with \(s^2\) \[Y \sim N(X\beta, \sigma^2 I)\] \[s^2 = \frac{e^{\prime}e}{n - p} = \frac{(Y - X\beta)^{\prime}(Y - X\beta)}{n - p}\] \[s^2 = \frac{SSE}{df_E} = MSE\]

• Variance-Covariance of Residuals

\[\epsilon = Y - \hat{Y} = Y - HY = (I - H) Y\]

\[Cov(\epsilon) = (I - H) Cov(Y) (I - H)^{\prime}\]

res <- y - X %*% beta

cbind(res, m$residuals) %>%
    prt(caption = "Residuals", col.names = (c("Matrix Algebra", "lm")))

sse <- (t(res) %*% res) / (n - p)

cbind(sqrt(sse), sigma(m)) %>%
    prt(caption = "Residual Mean Squared Errors (Sigma)", col.names = c("Matrix Algebra", "lm"))

Estimate Variance-Covariance of Coefficients

• Ref: https://timeseriesreasoning.com/contents/deep-dive-into-variance-covariance-matrices/
• Variance-covariance matrix of coefficients \[Var(\hat{\beta}) = \sigma^2 (X^{\prime}X)^{-1} = s^2 (X^{\prime}X)^{-1}\] \[V[\hat{\beta}] = E[\hat{\beta}^2] - E[\hat{\beta}^{\prime}] E[\hat{\beta}]\]

\[V[\hat{\beta}] = V[(X^TX)^{-1}X^Ty]\]

\[V[\hat{\beta}] = V[(X^TX)^{-1}X^T (X\hat{\beta} + \epsilon)] \] \[V[\hat{\beta}] = V[(X^TX)^{-1}X^T \epsilon]\] \[A = (X^TX)^{-1}X^{T}\] \[V[\hat{\beta}] = A V[\epsilon] A^T\]

A <- solve(t(X) %*% X) %*% t(X)
myCov <- A %*% diag(x = as.numeric(sse), nrow = n) %*% t(A)

cbind(myCov, vcov(m)) %>% prt(caption = "Variance-Covariance Matrix of Coefficients")

t-test on a single coeffient

• Test statistic \[t^* = (\hat{\beta} - 0)/ s(\beta)\] with degree of freedom \(n-p\)

t_values <- (beta - 0) / sqrt(diag(myCov))
p_values <- pt(t_values, df = n - p, lower.tail = FALSE) * 2

ses <- sqrt(diag(myCov))

cbind(Beta = beta, SE = ses, t_ = t_values, p_ = p_values, ms$coefficients) %>%
    prt(caption = "t-test on Coefficients")

Wald Statistic

• Ref: https://www.statisticshowto.com/wp-content/uploads/2016/09/2101f12WaldWithR.pdf
• Ref: https://cran.r-project.org/web/packages/clubSandwich/vignettes/Wald-tests-in-clubSandwich.html
• Ref: https://www.stat.umn.edu/geyer/8112/notes/tests.pdf
• Ref: http://home.cc.umanitoba.ca/~godwinrt/4042/material/part12.pdf
• Ref: https://bookdown.org/mike/data_analysis/wald-test.html
• This is Wald test statistic in a simple form \[W = (\hat{\theta} - \theta_0)^{\prime}[cov(\hat{\theta})]^{-1}(\hat{\theta} - \theta_0)\] \[W \sim \chi^2_q\] where \(q = rank(cov(\hat{\theta}))\)
• A linear restriction uses a set of values that make possible restrictions on \(\beta\) values
• \[R \beta = q\] is a restriction with dimensions \((j \times k)(k \times 1) = (j \times 1)\)
• If \(R = |1, 0, 0|\), then the restriction is that the first \(\beta_1 == q\)
• If \(R = |0, 1, 1|\), then the restriction is that the \(\beta_2 + \beta_3 == q\)
• Set \[m = R\beta - q\], then the null hypothesis is \(E[m] = 0\)
• And, \(V[m] = V[R\beta - q] = V[R\beta] = R V[\beta] R^{\prime}\)
• So, \(m \sim N[0, RV[\beta]R^{\prime}]\)
• The Wald test statistic \[W = m^{\prime} V[m]^{-1} m\]

R <- matrix(c(0, 1, 0, 0, 0, 1), nrow = 2, byrow = TRUE)
q <- matrix(c(0, 0), ncol = 1)
ctr <- R %*% beta - q
Vm <- R %*% myCov %*% t(R)
W <- t(ctr) %*% solve(Vm) %*% ctr

pchisq(q = W, df = dim(q)[1], lower.tail = FALSE)

     [,1]

[1,] 0.003458434

library(aod)
print(aod::wald.test(Sigma = myCov
                   , b = beta
                   , Terms = c(2:3)
                   , verbose = TRUE)
    , digits = 5)

Wald test:

Coefficients: [,1]
x0 2.81270 x1 1.53920 x2 0.33233

Var-cov matrix of the coefficients: x0 x1 x2
x0 2.121100 -0.544783 -0.047903 x1 -0.544783 0.337439 -0.037801 x2 -0.047903 -0.037801 0.063307

Test-design matrix: [,1] [,2] [,3] L1 0 1 0 L2 0 0 1

Positions of tested coefficients in the vector of coefficients: 2, 3

H0: 1.5392029 = 0; 0.3323287 = 0

Chi-squared test: X2 = 11.334, df = 2, P(> X2) = 0.0034584

95% Confidence Intervals for Coefficients

• \[CI = \hat{\beta} \pm t^c s(\beta)\]
• where \[t^c = t_{n-p}(0.975)\]

ci_lower <- beta - qt(0.975, df = n - p) * ses
ci_upper <- beta + qt(0.975, df = n - p) * ses
cbind(ci_lower, ci_upper, confint(m))

                      2.5 %    97.5 %

x0 -0.1755833 5.8009871 -0.1755833 5.8009871 x1 0.3473048 2.7311011 0.3473048 2.7311011 x2 -0.1839315 0.8485889 -0.1839315 0.8485889

Predicted Value

• Predict a subpopulation mean corresponding to a set of explanatory variables \[E(Y_i) = \mu_i = X_i^{\prime}\beta\] \[\hat{\mu_i} = X_i^{\prime} \hat{\beta}\] \[s^2(\hat{\mu}) = X^{\prime} s^2 (\beta) X = s^2 X^{\prime}(X^{\prime}X)^{-1}X\]
• 95% CI: \[\hat{\mu}_i \pm s(\hat{\mu}_i) t_{n-p}(0.975)\]

muhat <- X %*% beta



pred <- predict(m, newdata = dt, se.fit = TRUE, interval = "confidence")

semu <- diag(sqrt(diag(sse[1, 1, drop = TRUE], nrow = n) %*% X %*% solve(t(X) %*% X) %*% t(X)))

cbind(semu, pred$se.fit) %>% prt()

Predicted Mean

• Ref: http://www.math.kent.edu/~reichel/courses/monte.carlo/alt4.7d.pdf

\[E[\hat{y}] = E[X\hat{\beta}] = \sum (X\hat{\beta})/N\] \[V[E(\hat{y})] = V[\sum (X\hat{\beta})/N]\] \[V[E(\hat{y})] = \frac{1}{N^2} \sum(V[X\hat{\beta}])\] \[V[E(\hat{y})] = \frac{1}{N^2} \sum(XV[\hat{\beta}]X^T)\]

wt <- matrix(rep(1 / n, n), nrow = 1)
beta <- matrix(beta, ncol = 1)

mpred <- wt %*% (X %*% beta)

vpm <- diag(X %*% myCov %*% t(X))
wvpm <- wt %*% vpm

sqrt(sum(diag(vpm)))

[1] 9.677951

library(emmeans)
emm <- emmeans(object = m, specs = ~ 1)

print(emm, digits = 5)

1 emmean SE df lower.CL upper.CL overall 6.23 1.02 27 4.14 8.32

Confidence level used: 0.95

Weighted Predicted Mean

\[E[w\hat{y}] = E[wX\beta] = \sum (wX\beta)/N\] \[V[w\hat{y}] = V[wX\beta] = V[\sum (wX\beta)/N]\] \[V[w\hat{y}] = \frac{1}{N^2} \sum(V[wX\beta])\] \[V[w\hat{y}] = \frac{1}{N^2} \sum(w^2V[X\beta])\]

wt <- matrix(dt$w, nrow = 1)
beta <- matrix(beta, ncol = 1)

wpred <- wt %*% (X %*% beta)

Robust Standard Error

• The errors terms have constant variance and are uncorrelated \[E[uu^T] = \sigma^2 I_n\]
• When variance is not constant \[V[\hat{\beta}_{OLS}] = V[(X^TX)^{-1}X^Ty] = (X^TX)^{-1}X^T \sum X(X^TX)^{-1}\] where \(\sum = V[u]\)
• White’s method, or HCE (heteroskedasticity-consistent estimator), or HC0 \[\hat{V}_{HCE}[\hat{\beta}_{OLS}] = (X^TX)^{-1}(X^T diag(\hat{\epsilon}_1^2, ..., \hat{\epsilon}_n^2) X)(X^TX)^{-1}\]

sgm <- res ^ 2
sgm <- diag(as.vector(sgm), nrow = dim(sgm)[1])

vhce <- solve(t(X) %*% X) %*% (t(X) %*% sgm %*% X) %*% solve(t(X) %*% X)



library(sandwich)
rbind(vhce, vcovHC(m, type = "HC0")) %>% prt()

Model Matrix

set.seed(123)
mydata <- data.frame(
  response = rnorm(10, mean = 50, sd = 10),
  predictor = rnorm(10, mean = 50, sd = 10),
  group = as.factor(rep(1:2, each = 5))
)

fm1 <- lm(response ~ predictor + group, data = mydata)
fm2 <- lme4::lmer(response ~ predictor + (1|group), data = mydata)
fm3 <- nlme::lme(response ~ predictor, random = ~1|group, data = mydata)

summary(fm1)

Call: lm(formula = response ~ predictor + group, data = mydata)

Residuals: Min 1Q Median 3Q Max -12.3400 -3.7151 -0.5807 3.6658 13.1649

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 24.1296 15.5672 1.550 0.165 predictor 0.5239 0.2838 1.846 0.107 group2 -1.3387 5.5892 -0.240 0.818

Residual standard error: 8.792 on 7 degrees of freedom Multiple R-squared: 0.3391, Adjusted R-squared: 0.1502 F-statistic: 1.795 on 2 and 7 DF, p-value: 0.2347

summary(fm2)

Linear mixed model fit by REML [‘lmerMod’] Formula: response ~ predictor + (1 | group) Data: mydata

REML criterion at convergence: 65.7

Scaled residuals: Min 1Q Median 3Q Max -1.42166 -0.40578 -0.06708 0.48858 1.67363

Random effects: Groups Name Variance Std.Dev. group (Intercept) 0.0 0.000
Residual 68.2 8.258
Number of obs: 10, groups: group, 2

Fixed effects: Estimate Std. Error t value (Intercept) 23.1034 14.0566 1.644 predictor 0.5307 0.2652 2.001

Correlation of Fixed Effects: (Intr) predictor -0.983 optimizer (nloptwrap) convergence code: 0 (OK) boundary (singular) fit: see help(‘isSingular’)

dim(model.matrix(fm1))

[1] 10 3

model.matrix(fm1)[1:3, ]

(Intercept) predictor group2 1 1 62.24082 0 2 1 53.59814 0 3 1 54.00771 0

getME(fm2, "X")[, ]

(Intercept) predictor 1 1 62.24082 2 1 53.59814 3 1 54.00771 4 1 51.10683 5 1 44.44159 6 1 67.86913 7 1 54.97850 8 1 30.33383 9 1 57.01356 10 1 45.27209

getME(fm2, "Z")[, ]

10 x 2 sparse Matrix of class “dgCMatrix” 1 2 1 1 . 2 1 . 3 1 . 4 1 . 5 1 . 6 . 1 7 . 1 8 . 1 9 . 1 10 . 1

model.matrix(terms(fm3), data = mydata)

(Intercept) predictor 1 1 62.24082 2 1 53.59814 3 1 54.00771 4 1 51.10683 5 1 44.44159 6 1 67.86913 7 1 54.97850 8 1 30.33383 9 1 57.01356 10 1 45.27209 attr(,“assign”) [1] 0 1

model.matrix(terms(fm3, random = TRUE), data = mydata)

(Intercept) predictor 1 1 62.24082 2 1 53.59814 3 1 54.00771 4 1 51.10683 5 1 44.44159 6 1 67.86913 7 1 54.97850 8 1 30.33383 9 1 57.01356 10 1 45.27209 attr(,“assign”) [1] 0 1

model.matrix(attr(fm3$modelStruct$reStr$group, "formula"), data = mydata)

(Intercept) 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 1 10 1 attr(,“assign”) [1] 0

model.matrix(formula(fm3$modelStruct$reStr)[[1]], data = mydata)

(Intercept) 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 1 10 1 attr(,“assign”) [1] 0

Compound-Symmetry Structure

• In an example to model slope over multiple times on same individual
• Introducing a common correlation among the observations from a single individual

• In additiona to \(\sigma^2\), one CS covariance \(\sigma_1^2\) is added

• CS sturcture in R matrix without setting Z and G

• This is the covariance structure for one person assuming three repeated measurements

\[\left[\begin{matrix} \sigma_1^2 + \sigma^2 & \sigma_1^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 + \sigma^2 & \sigma_1^2 \\ \sigma_1^2 & \sigma_1^2 & \sigma_1^2 + \sigma^2 \end{matrix}\right]\]

• This is coded in SAS as

proc mixed;
class indiv;
model y = time;
repeated / type=cs subject=indiv;
run;

• This could also be coded in SAS with G and Z matrix

proc mixed;
class indiv;
model y = time;
random indiv;
run;

• Another way of coding in SAS with G and Z matrix

proc mixed;
class indiv;
model y = time;
random intercept / subject=indiv;
run;

Two Random Effects

• This is SAS coding

proc mixed;
class a block;
model y = a;
random block a*block;
run;

• Alternative SAS coding

random int a / subject=block;

• Assume a has two levels and block has 3 levels, then each individual has 6 repeated measurements
• This is an example of random effect on the block level
• Z matrix, the first two columns are of a levels, and the last three columns are of block levels

\[\left[\begin{matrix} 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 1 \end{matrix}\right]\]

• G matrix where \(\sigma_B^2\) is the variance component for block, this is the default VC (variance components) for random statement in SAS by adding one variance for each random effect

\[\left[\begin{matrix} \sigma_B^2 & & & & & \\ & \sigma_B^2 & & & & \\ & & \sigma_B^2 & & & \\ & & & \sigma_{AB}^2 & & \\ & & & & \sigma_{AB}^2 & \\ & & & & & \sigma_{AB}^2 \end{matrix}\right]\]