Statistics 892 Linear Algebra/Models Readings - Assignment 1

1. View videos Stat902-012317; Stat902-012517; Stat902-013017 from box.unl.edu: https://unl.app.box.com/folder/66083066660?s=idup7tp7rvkh8w3fxjqpo52u6f2xzc9m

2. Read Johnson and Wichern (6th ed) JW, Chapters 2.1-2.2; 2.5-2.6 and read Ch. 2 supplement as needed. Work exercises 2.2, 2.5, 2.7(c) only.

2.2

A = matrix(c(-1,4,3,2), nrow = 2)
B = matrix(c(4,1,-2,-3,-2,0), nrow = 3)
C = matrix(c(5,-4,2), ncol=1)

#5A
A5 = A*5
A5
##      [,1] [,2]
## [1,]   -5   15
## [2,]   20   10
#BA
BA = B%*%A
BA  
##      [,1] [,2]
## [1,]  -16    6
## [2,]   -9   -1
## [3,]    2   -6
#A'B'
Ap = t(A)
Bp = t(B)
ApBp = Ap%*%Bp
ApBp
##      [,1] [,2] [,3]
## [1,]  -16   -9    2
## [2,]    6   -1   -6
#C'B
Cp = t(C)
CpB = Cp%*%B
CpB
##      [,1] [,2]
## [1,]   12   -7
#is AB defined?
#no it is a non-conformable argument

2.7(c)

A = matrix(c(9,-2,-2,6), nrow = 2)
A.inv = solve(A)
A
##      [,1] [,2]
## [1,]    9   -2
## [2,]   -2    6
q = matrix(c(5/13, -12/13, 12/13, 5/13), nrow = 2)
q.transpose = t(q)
I = q.transpose%*%q
I
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
#the identity matrix is produced and indicates the matrix is orthogonal

3. Read JW Ch7.1-7.3 thru p367; 7.4 and example thru 374. Work exercise 7.1 using matrix operations in R or SAS IML and in addition

Z = matrix(c(1,10,1,5,1,7,1,19,1,11,1,8), nrow = 2)
y = matrix(c(15, 9, 3, 25, 7, 13), nrow =1)
# transpose of z = z'
Zp = t(Z) 
# beta.hat = (z'z)^-1*z'y
beta= y%*%(Zp%*%(solve(Z%*%Zp))) 
beta #least squares solution vector
##            [,1]     [,2]
## [1,] -0.6666667 1.266667
#vector of fitted predicted values
y.hat = beta%*%Z
y.hat
##      [,1]     [,2] [,3] [,4]     [,5]     [,6]
## [1,]   12 5.666667  8.2 23.4 13.26667 9.466667
#residual sum of squares ê'e
e = y - y.hat
epe = e%*%t(e)
epe
##          [,1]
## [1,] 101.4667

(b) compute the standard error of the estimated slope coefficient and interpret its meaning,

Z = matrix(c(10,5,7,19,11,8), nrow = 1)
Y = matrix(c(15,9,3,25,7,13), nrow = 1)
SZY <- sum((Z-mean(Z))*(Y-mean(Y)))
SZZ <- sum((Z-mean(Z))^2)
SYY <- sum((Y-mean(Y))^2)
beta.hat.1 <- SZY/SZZ
beta.hat.1
## [1] 1.266667
beta.hat.0 <- mean(Y) - beta.hat.1 * mean(Z)
beta.hat.0
## [1] -0.6666667
e.hat <- Y - (beta.hat.0 + beta.hat.1 * Z)
e.hat
##      [,1]     [,2] [,3] [,4]      [,5]     [,6]
## [1,]    3 3.333333 -5.2  1.6 -6.266667 3.533333
RSS <- sum(e.hat^2)
RSS
## [1] 101.4667
n <- length(Y)
sigma.hat.squared <- RSS / (n-2)
sigma.hat.squared
## [1] 25.36667
sigma.hat <- sqrt(sigma.hat.squared)
sigma.hat
## [1] 5.036533

(c) test the hypothesis that there is no relationship between z1 and y.

Ho: µ1 = µ2, where there is no relationship between z1 and y

Ha: µ1 ≠ µ2, where there is a relationship between z1 and y

observations <- c(10,5,7,19,11,8,15,9,3,25,7,13)
treatments <- c("z","z","z","z","z","z","y","y","y","y","y","y")
data <- cbind(observations,treatments)
colnames(data) <- c("obs", "trt")
data <- data.frame(data)
data$obs <- as.numeric(data$obs)
data$trt <- as.factor(data$trt)

model <- aov(obs~trt, data = data)
summary(model)
##             Df Sum Sq Mean Sq F value Pr(>F)
## trt          1     12    12.0    0.29  0.602
## Residuals   10    414    41.4
summary.lm(model)
## 
## Call:
## aov(formula = obs ~ trt, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -9.0   -3.5   -1.0    1.5   13.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   12.000      2.627   4.568  0.00103 **
## trtz          -2.000      3.715  -0.538  0.60209   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.434 on 10 degrees of freedom
## Multiple R-squared:  0.02817,    Adjusted R-squared:  -0.06901 
## F-statistic: 0.2899 on 1 and 10 DF,  p-value: 0.6021

the p-value = 0.602 shows that we fail to reject the Ho and concur that there is no relationship between z1 and y

(d) compute the model and error SS (SSmodel and SSE) and

I can answer this by computing through the anova() statement

model1 <- anova(model)
SSmodel <- model1$'Sum Sq'[1] + model1$'Sum Sq'[2] 
SSmodel
## [1] 426
SSE <- sum(model$residuals^2)
SSE
## [1] 414

Or I can go through the mathematical steps

Z <- matrix(c(10,5,7,19,11,8), nrow = 1)
Y <- matrix(c(15,9,3,25,7,13), nrow = 1)
#total number of observations
n <- length(Z)+length(Y)
#length of Z and Y
Zn <- length(Z)
Yn <- length(Y)
#number of treatments - Z, Y
k <- 2
#mean of Z and Y
Zµ <- mean(Z)
Zµ
## [1] 10
Yµ <- mean(Y)
Yµ
## [1] 12
#grand mean
Gµ <- mean(c(Zµ, Yµ))
Gµ
## [1] 11
#sum of squares ZY
SZY <- sum((Z-mean(Z))*(Y-mean(Y)))
#sum of squares ZZ
SZZ <- sum((Z-mean(Z))^2)
#sum of squares YY
SYY <- sum((Y-mean(Y))^2)
#sum of squares within
SSw <- SZZ + SYY
SSw #Residuals Sum Sq
## [1] 414
#mean sum of squares within
MSSw <- SSw/(n-k)
MSSw #Residuals Mean Sq
## [1] 41.4
#sum of squares between
SSb <- sum((Zn*((Zµ - Gµ)^2)) + (Yn*((Yµ - Gµ)^2)))
SSb #trt Sum Sq
## [1] 12
#mean sum of squares between
MSSb <- SSb/(k-1)
MSSb #obs Sum Sq 
## [1] 12
#compute F-statistic
Fstat = MSSb/MSSw
Fstat #F value
## [1] 0.2898551
#compute p-value from Fstat that uses df from trt (df = 2-1) and Residuals (df = 12-2)
pf(Fstat, 1, 10, lower.tail = FALSE)
## [1] 0.6020904

(e) compute R2 as 1 – Ssmodel / SSE and interpret it’s meaning.

Using R functions I can calculate the coefficient of determination, R^2

R <- cor(observations, fitted(model))
R.squared <- R^2
R.squared
## [1] 0.02816901

Or, I can calculate manually using SSw and SSb, where SSw is the residual

R.squared = 1 - SSw/(SSw+SSb)
R.squared
## [1] 0.02816901

4. An experiment was conducted as a completely randomized design (CRD) with a single treatment factor with three levels. Two experimental units were used per treatment level. The following responses were obtained:

TRT —-1—- —–2—- —-3—-

Y 10 12 14 16 9 8

a. Use R lm() or aov() OR SAS GLM or MIXED or GLIMMIX (or another statistical program) to obtain the obtain the least squares solution vector. In SAS, use the SOLUTION option in the MODEL statement.

y = matrix(c(10,12,14,16,9,8),nrow=6)
trt = c(1,1,2,2,3,3)
df = data.frame(cbind(y,trt))
colnames(df)[1] <- "score"
df$trt = factor(df$trt)
options(contrasts = c("contr.treatment", "contr.poly"))
fit = aov(score ~ trt, data = df)
summary(fit)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## trt          2   43.0    21.5   14.33 0.0292 *
## Residuals    3    4.5     1.5                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.lm(fit)
## 
## Call:
## aov(formula = score ~ trt, data = df)
## 
## Residuals:
##    1    2    3    4    5    6 
## -1.0  1.0 -1.0  1.0  0.5 -0.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   11.000      0.866  12.702  0.00105 **
## trt2           4.000      1.225   3.266  0.04692 * 
## trt3          -2.500      1.225  -2.041  0.13389   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.225 on 3 degrees of freedom
## Multiple R-squared:  0.9053, Adjusted R-squared:  0.8421 
## F-statistic: 14.33 on 2 and 3 DF,  p-value: 0.02916
x = model.matrix(fit)
x
##   (Intercept) trt2 trt3
## 1           1    0    0
## 2           1    0    0
## 3           1    1    0
## 4           1    1    0
## 5           1    0    1
## 6           1    0    1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$trt
## [1] "contr.treatment"
beta = solve(t(x)%*%x)%*%t(x)%*%y
beta # least squares solution vector
##             [,1]
## (Intercept) 11.0
## trt2         4.0
## trt3        -2.5

Also, test the hypothesis (μ1+μ2)/2 = μ3.. In SAS use the CONTRAST statement. If you are using R, with lm() or aov(), make sure to identify TRT as a factor using the factor() function and use either the fit.contrast() or glht() functions.) To compute parts b and c, use matrix calculationswith either R or SAS PROC IML, or another package (eg Matlab, Python).

b. Set up the design matrix (“X”) for the problem using the “set to zero” restriction with τ3= 0.

x = t(matrix(c(1,1,0,1,1,0,1,0,1,1,0,1,1,0,0,1,0,0), nrow = 3))
x
##      [,1] [,2] [,3]
## [1,]    1    1    0
## [2,]    1    1    0
## [3,]    1    0    1
## [4,]    1    0    1
## [5,]    1    0    0
## [6,]    1    0    0

c. Use the matrix statements in R or PROC IML (or another statistical package) to obtain the following. Cut and paste the relevant portion of the output into your homework paper. In SAS, your solution should be the same from the SAS solution. It will not be equal to the lm() or aov() solution vector. Explain why SAS and R differ regarding the solutions.

the main difference between SAS and R is how they set a t = 0, in R it sets the last t to zero, where sas sets the first to zero.

(1) Compute the least squares solution using the τ3= 0 set to zero restriction.

to calculate beta, this is using the x matrix created in (b).

beta = solve(t(x)%*%x)%*%t(x)%*%y
beta # least squares solution vector
##      [,1]
## [1,]  8.5
## [2,]  2.5
## [3,]  6.5

(2) Estimate σ^2.

#using the model to extract sigma and square it
sigma <- sigma(fit, use.fallback = TRUE)
sigma # residual standard error
## [1] 1.224745
sigma.hat.squared1 <- sigma^2
sigma.hat.squared1
## [1] 1.5
#using matrix statement to calculate sigma.hat.squared, x matrix of t3 = 0
sigma.hat.squared <- as.numeric((t(y)%*%y - t(beta)%*%t(x)%*%y)/(length(y) - length(unique(trt))))
sigma.hat.squared
## [1] 1.5

(3) Estimate the covariance matrix of the least squares solution vector obtained from part (1) above.

cov.beta.hat <- sigma.hat.squared*(solve(t(x)%*%x))
cov.beta.hat
##       [,1]  [,2]  [,3]
## [1,]  0.75 -0.75 -0.75
## [2,] -0.75  1.50  0.75
## [3,] -0.75  0.75  1.50

(4) Using the covariance matrix from (3), compute the standard errors of the estimates.

The diagonal of the square root of the cov.beta.hat matrix

se.cov <- sqrt(diag(cov.beta.hat))
se.cov
## [1] 0.8660254 1.2247449 1.2247449

(5) Obtain the total corrected sums of squares, the sums of squares for treatment (i.e. R(τ1,τ2,τ3|μ) ) and the sums of squares for error.

#using model statements
fit1 <- anova(fit)
ss.total <- fit1$'Sum Sq'[1] + fit1$'Sum Sq'[2]
ss.total
## [1] 47.5
ss.trt <- fit1$'Sum Sq'[1]
ss.trt
## [1] 43
ss.error <- fit1$'Sum Sq'[2]
ss.error
## [1] 4.5
#manually calculating
#df1 <- matrix(c(10,12,14,16,9,8),ncol=3)
#df1 <- data.frame(df1)
#colnames(df1) <- c("trt1","trt2","trt3")

#ssw.1 <- sum((df1$trt1-mean(c(df1$trt1)))^2)
#ssw.2 <- sum((df1$trt2-mean(c(df1$trt2)))^2)
#ssw.3 <- sum((df1$trt3-mean(c(df1$trt3)))^2)
#grand.mean <- mean(c(mean(df$trt1), mean(df$trt2), mean(df$trt3)))
#SSw <- ssw.1+ssw.2+ssw.3 #sse
#SSw

#ssb.1 <- length(df1$trt1)*((mean(c(df$trt1, na.rm = TRUE)) - grand.mean)^2)
#ssb.2 <- length(df1$trt2)*((mean(c(df$trt2, na.rm = TRUE)) - grand.mean)^2)
#ssb.3 <- length(df1$trt3)*((mean(c(df$trt3, na.rm = TRUE)) - grand.mean)^2)
#SSb <- ssb.1+ssb.2+ssb.3
#SSb
#TSS <- SSw + SSb
#TSS

(6) Develop a K’ matrix and test H0:(μ1+μ2)/2 = μ3 as described in the first example in the Linear Models section of the first handout. Note that μi = μ + τi i=1,2,3.

xpx <- t(x)%*%x
e <- (y-x%*%beta)
e
##      [,1]
## [1,] -1.0
## [2,]  1.0
## [3,] -1.0
## [4,]  1.0
## [5,]  0.5
## [6,] -0.5
epe <- t(e)%*%e
epe
##      [,1]
## [1,]  4.5
mse <- epe/3
mse
##      [,1]
## [1,]  1.5
kp <- t(matrix(c(0,1,-1,0,1,0),ncol=2))
kp
##      [,1] [,2] [,3]
## [1,]    0    1   -1
## [2,]    0    1    0
r <- matrix(c(0,0))
f <- t(kp%*%beta - r)%*%solve(kp%*%solve(xpx)%*%t(kp))%*%(kp%*%beta - r) / (2%*%mse)
f
##          [,1]
## [1,] 14.33333
pvalue <- 1- pf(f,2,3)
pvalue
##            [,1]
## [1,] 0.02915938
kp <- t(matrix(c(1, -1/2, -1/2), ncol = 1))
kp
##      [,1] [,2] [,3]
## [1,]    1 -0.5 -0.5
r <- matrix(c(0))
f1 <- t(kp%*%beta - r)%*%solve(kp%*%solve(xpx)%*%t(kp))%*%(kp%*%beta - r) / (1%*%mse)
f1
##          [,1]
## [1,] 4.740741
pvalue <- 1 - pf(f,2,3)
pvalue
##            [,1]
## [1,] 0.02915938
#contrasts
ll <- aov(score ~ trt, data = df)
summary(ll);
##             Df Sum Sq Mean Sq F value Pr(>F)  
## trt          2   43.0    21.5   14.33 0.0292 *
## Residuals    3    4.5     1.5                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(gmodels)
trt1vs2and3 <- fit.contrast(ll, "trt", c(1, -1/2, -1/2)) #(µ1+µ2)/2 = µ3
trt1vs2and3
##                       Estimate Std. Error    t value  Pr(>|t|)
## trt c=( 1 -0.5 -0.5 )    -0.75    1.06066 -0.7071068 0.5304778
## attr(,"class")
## [1] "fit_contrast"
trt1vs2 <- fit.contrast(ll, "trt", c(1, -1, 0))
trt1vs2
##                  Estimate Std. Error   t value   Pr(>|t|)
## trt c=( 1 -1 0 )       -4   1.224745 -3.265986 0.04691893
## attr(,"class")
## [1] "fit_contrast"
trt1vs3 <- fit.contrast(ll, "trt",c(1, 0, -1))
trt1vs3
##                  Estimate Std. Error  t value  Pr(>|t|)
## trt c=( 1 0 -1 )      2.5   1.224745 2.041241 0.1338866
## attr(,"class")
## [1] "fit_contrast"
#using multiple comparisons of means test
library(multcomp);
k <- rbind(c(1, -1/2, -1/2),
           c(1, -1, 0),
           c(1, 0, -1))
fit.gh <- glht(ll, linfct = mcp(trt = k))
summary(fit.gh, test = adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = score ~ trt, data = df)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)  
## 1 == 0   -0.750      1.061  -0.707   0.5305  
## 2 == 0   -4.000      1.225  -3.266   0.0469 *
## 3 == 0    2.500      1.225   2.041   0.1339  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)