forked from vasishth/MScStatisticsNotes
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathExampleslinmod.Rnw
90 lines (65 loc) · 1.36 KB
/
Exampleslinmod.Rnw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
\documentclass{article}
\begin{document}
\SweaveOpts{concordance=TRUE}
\title{Linear modelling critical examples}
\section{Hypothesis testing using C matrix}
<<>>=
data(ToothGrowth)
fm1<-lm(len~dose+supp,ToothGrowth)
#summary(fm1)
## see below:
(varcov.mat<-round(vcov(fm1),digits=3))
@
Test: $H_0: \beta_0=2\beta_1, \beta_2=0$.
<<>>=
y<-ToothGrowth$len
X<-model.matrix(fm1)
XT.X <- t (X)%*%X
XT.y <- t(X)%*%y
G<-solve(XT.X)
## recovers Vcov matrix:
4.24^2*G
beta.hat <- G%*%XT.y
C<-matrix(c(1,-2,0,0,0,1),byrow=T,ncol=3)
## this changes:
c<-matrix(c(0,0),byrow=T,ncol=1)
(yT.y <- t(y)%*%y)
(bT.XT.y<-t(beta.hat)%*%t(X)%*%y)
## by definition of S_r:
Sr <- yT.y - bT.XT.y
n<-60
p<-3 ## num parameters
sigma.hat.2 <- (1/(n-p)) * Sr
#sqrt(sigma.hat.2)
num<-t(C %*% beta.hat - c) %*%
ginv(C %*% G %*% t(C) ) %*%
(C %*% beta.hat - c )
q<-2 ## num of hypotheses
(F<-num/(sigma.hat.2*q))
1-pf(F,q,n-p)
@
In the exam exercise we are given:
<<>>=
given<-solve(C%*%solve(XT.X)%*%t(C))/4.24^2
@
So, the middle term is:
<<>>=
given*4.24^2
@
Cross-check:
<<>>=
ginv(C %*% G %*% t(C) )
@
So we can do:
<<>>=
num<-t(C %*% beta.hat - c) %*%
ginv(C %*% G %*% t(C) ) %*%
(C %*% beta.hat - c )
num2<-t(C %*% beta.hat - c) %*%
(given*4.24^2)%*%
(C %*% beta.hat - c )
q<-2 ## num of hypotheses
(F<-num/(sigma.hat.2*q))
1-pf(F,q,n-p)
@
\end{document}