回歸作業

Y=c(78.5,74.3,104.3,87.6,95.9,109.2,102.7)
X1=c(7,1,11,11,7,11,3)
X2=c(2.6,2.9,5.6,3.1,5.2,5.5,7.1)
X0=c(rep(1,7))

#a
X=cbind(X0,X1,X2)
X

##      X0 X1  X2
## [1,]  1  7 2.6
## [2,]  1  1 2.9
## [3,]  1 11 5.6
## [4,]  1 11 3.1
## [5,]  1  7 5.2
## [6,]  1 11 5.5
## [7,]  1  3 7.1

#b
t(X)%*%X

##    X0  X1     X2
## X0  7  51  32.00
## X1 51 471 235.00
## X2 32 235 163.84

t(X)%*%Y

##       [,1]
## X0  652.50
## X1 4915.30
## X2 3103.66

#c
solve(t(X)%*%X)

##             X0          X1          X2
## X0  1.79959717 -0.06854721 -0.25316476
## X1 -0.06854721  0.01007738 -0.00106613
## X2 -0.25316476 -0.00106613  0.05707894

beta.hat<-solve(t(X)%*%X)%*%t(X)%*%Y
model<-lm(Y~X1+X2)
model$coefficients

## (Intercept)          X1          X2 
##   51.569698    1.497410    6.723256

#Y=51.569698+1.497410*X1+6.723256*X2
#當X1每增加1時，Y會增加1.497410
#當X2每增加1時，Y會增加6.723256

#d
Y.hat<-X%*%beta.hat
e=Y-Y.hat
Y.hat

##           [,1]
## [1,]  79.53204
## [2,]  72.56455
## [3,] 105.69144
## [4,]  88.88331
## [5,]  97.01250
## [6,] 105.01912
## [7,] 103.79704

e

##           [,1]
## [1,] -1.032036
## [2,]  1.735450
## [3,] -1.391445
## [4,] -1.283305
## [5,] -1.112501
## [6,]  4.180881
## [7,] -1.097045

#e
MSE=2.626
P=X%*%solve(t(X)%*%X)%*%t(X)
var.beta.hat=solve(t(X)%*%X)*MSE^2
var.y.hat=P*MSE^2
var.e=(diag(1,7)-P)*MSE^2

#第二題
library(alr3)

## Loading required package: car

x=c(0,0,0,4,4,4,6,6,6,6,6,8,8,8,8,12,12)
y=c(121,130,100,128,124,132,141,145,150,148,151,152,146,149,133,125,121)
df=data.frame(x,y)
kk=lm(y~x)
pureErrorAnova(kk)

## Analysis of Variance Table
## 
## Response: y
##              Df  Sum Sq Mean Sq F value   Pr(>F)   
## x             1  441.16  441.16  6.7012 0.023716 * 
## Residuals    15 2875.78  191.72                    
##  Lack of fit  3 2085.78  695.26 10.5609 0.001102 **
##  Pure Error  12  790.00   65.83                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1