LMT: Weeks 11 and 12 Examples
Heather Kitada Smalley
Delivery Data
delData<-read.csv("https://raw.githubusercontent.com/kitadasmalley/sp21_MATH376LMT/main/data/deliveryTime.csv",
header=TRUE)
attach(delData)Class 11A:
Create a 3D Scatterplot
#install.packages("scatterplot3d")
library("scatterplot3d")
sp3d<-scatterplot3d(x=numCases, y=distance, z=delTime,
pch=16, angle = 60)
Add the MLR regression plane
# Add regression plane
mod <- lm(delTime~numCases+distance)
#sp3d$plane3d(mod)Produce the H-matrix (projection matrix)
# Design mat
xMat<-matrix(c(rep(1, length(delTime)),
numCases, distance),
nrow=length(delTime))
head(xMat)## [,1] [,2] [,3]
## [1,] 1 7 560
## [2,] 1 3 220
## [3,] 1 3 340
## [4,] 1 4 80
## [5,] 1 6 150
## [6,] 1 7 330# Hat
hMat<-xMat%*%solve(t(xMat)%*%xMat)%*%t(xMat)Confirm that H(I-H)=0
# Identity matrix
I<-diag(length(delTime))
head(hMat%*%(I-hMat))## [,1] [,2] [,3] [,4] [,5]
## [1,] 8.673617e-18 4.987330e-18 -1.149254e-17 1.821460e-17 1.257675e-17
## [2,] -8.803722e-17 -1.474515e-17 -4.857226e-17 2.168404e-17 2.602085e-17
## [3,] -4.119968e-17 3.469447e-18 -2.168404e-17 2.038300e-17 1.474515e-17
## [4,] -1.088539e-16 -3.209238e-17 -8.196568e-17 1.561251e-17 3.295975e-17
## [5,] -9.324139e-17 -4.683753e-17 -8.196568e-17 3.556183e-17 2.688821e-17
## [6,] -4.770490e-17 -8.239937e-18 -4.076600e-17 1.951564e-17 1.908196e-17
## [,6] [,7] [,8] [,9] [,10]
## [1,] 8.890458e-18 1.474515e-17 1.409463e-17 -5.117434e-17 1.073360e-17
## [2,] -1.561251e-17 1.994932e-17 1.344411e-17 -1.578598e-16 -1.103718e-16
## [3,] -1.301043e-17 5.204170e-18 1.301043e-17 -1.860491e-16 -7.394259e-17
## [4,] -2.905662e-17 4.336809e-18 2.515349e-17 -1.812786e-16 -1.543904e-16
## [5,] -1.257675e-17 -8.673617e-18 3.252607e-17 -1.127570e-16 -1.262011e-16
## [6,] -5.637851e-18 1.387779e-17 1.821460e-17 -6.982262e-17 -8.239937e-17
## [,11] [,12] [,13] [,14] [,15]
## [1,] 3.198396e-18 9.324139e-18 8.239937e-18 8.456777e-18 2.168404e-18
## [2,] -5.063224e-17 3.382711e-17 -1.517883e-17 -6.028164e-17 -2.298509e-17
## [3,] -4.141652e-17 1.517883e-17 4.336809e-18 -3.642919e-17 -2.840610e-17
## [4,] -7.329207e-17 5.724587e-17 -2.862294e-17 -8.586881e-17 -4.900594e-17
## [5,] -1.517883e-17 6.982262e-17 -2.125036e-17 -7.589415e-17 -2.385245e-17
## [6,] -1.572093e-17 4.336809e-17 -1.517883e-17 -3.534499e-17 -1.105886e-17
## [,16] [,17] [,18] [,19] [,20]
## [1,] -1.593777e-17 1.019150e-17 2.038300e-17 1.778092e-17 -1.441989e-17
## [2,] -1.232196e-16 1.474515e-17 4.250073e-17 4.770490e-17 -7.214010e-17
## [3,] -7.583994e-17 1.301043e-18 2.081668e-17 2.341877e-17 -3.962759e-17
## [4,] -1.639314e-16 2.602085e-17 4.336809e-17 3.122502e-17 -5.076777e-17
## [5,] -1.407837e-16 3.469447e-18 5.551115e-17 1.127570e-17 -5.971243e-17
## [6,] -6.749159e-17 1.517883e-17 3.382711e-17 2.688821e-17 -2.469270e-17
## [,21] [,22] [,23] [,24] [,25]
## [1,] 1.105886e-17 1.119439e-17 1.517883e-18 2.331035e-17 6.938894e-18
## [2,] 5.811324e-17 -2.846031e-17 -4.553649e-17 -1.012645e-16 1.387779e-17
## [3,] 3.035766e-17 -1.246832e-17 -2.276825e-17 -5.225854e-17 6.938894e-18
## [4,] 9.280771e-17 3.480289e-17 -5.811324e-17 -1.292369e-16 1.387779e-17
## [5,] 8.760354e-17 2.569559e-17 -2.732189e-17 -9.952976e-17 0.000000e+00
## [6,] 6.028164e-17 2.688821e-17 -1.214306e-17 -5.724587e-17 6.938894e-18head(round(hMat%*%(diag(length(delTime))-hMat),2))## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0Class 11B:
Estimate the nuissance parameter \(\sigma\)
n<-length(delTime)
p<-2
ms_res<-sum(mod$residuals^2)/(n-p-1)
sigma_hat<-sqrt(ms_res)Estimate the covariance matrix for the \(\beta\) estimates
What are the standard errors for \(\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2\)?
What is the covariance between \(\hat{\beta}_1\) and \(\hat{\beta}_2\)
covMat<-ms_res*solve(t(xMat)%*%xMat)
sqrt(diag(covMat))## [1] 1.096730168 0.170734918 0.003613086# Confirm with
summary(mod)##
## Call:
## lm(formula = delTime ~ numCases + distance)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7880 -0.6629 0.4364 1.1566 7.4197
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.341231 1.096730 2.135 0.044170 *
## numCases 1.615907 0.170735 9.464 3.25e-09 ***
## distance 0.014385 0.003613 3.981 0.000631 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.259 on 22 degrees of freedom
## Multiple R-squared: 0.9596, Adjusted R-squared: 0.9559
## F-statistic: 261.2 on 2 and 22 DF, p-value: 4.687e-16Class 12A:
Use quadratic forms to define SS(Tot), SS(Reg), and SS(Res)
# Response vect
Y<-delTime
### SS(Res) = t(Y)%*%(diag(length(delTime))-hMat)%*%Y
I<-diag(length(delTime))
ss_res<-t(Y)%*%(I-hMat)%*%Y
ss_res## [,1]
## [1,] 233.7317### SS(Tot)
ones<-matrix(rep(1, length(delTime)),
nrow=length(delTime))
J<-ones%*%t(ones)
ss_tot<-t(Y)%*%(I-(1/n)*J)%*%Y
ss_tot## [,1]
## [1,] 5784.543### SS(Reg)
ss_reg<-t(Y)%*%(hMat-(1/n)*J)%*%Y
ss_reg## [,1]
## [1,] 5550.811Put together an ANOVA table
### MS(Res)
ms_res<-ss_res/(n-p-1)
ms_res## [,1]
## [1,] 10.62417### MS(Reg)
ms_reg<-ss_reg/p
ms_reg## [,1]
## [1,] 2775.405### F-stat
f_stat<-ms_reg/ms_res
f_stat## [,1]
## [1,] 261.2351### P-value
pf(f_stat, df1=p, df2=(n-p-1), lower.tail = FALSE)## [,1]
## [1,] 4.687422e-16# reject the null in favor of the alternative
## COMPARE THIS TO anova()
# note that there are seperate lines for the different explanatory variables
anova(mod)## Analysis of Variance Table
##
## Response: delTime
## Df Sum Sq Mean Sq F value Pr(>F)
## numCases 1 5382.4 5382.4 506.619 < 2.2e-16 ***
## distance 1 168.4 168.4 15.851 0.0006312 ***
## Residuals 22 233.7 10.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Look at cutt-offs from the F-distribution
### Calculate cuttoffs (we always want this upper tail for an F-test)
qf(0.01, df1=p, df2=(n-p-1), lower.tail = FALSE)## [1] 5.719022qf(0.05, df1=p, df2=(n-p-1), lower.tail = FALSE)## [1] 3.443357qf(0.1, df1=p, df2=(n-p-1), lower.tail = FALSE)## [1] 2.561314detach(delData)