Gumbel: Exercise 1

From MathWiki

Table of contents

Male:1, Female:0

Without Interaction Terms

> ### read input and make a copy ###
> dd <- read.table("data1.txt", header = T)
> dd2 <- dd
> 
> ### make required variables ###
> dd2$Sex <- ifelse(dd2$Sex == "M",1,0)
> dd2$sex.Educ <- dd2$Sex*dd2$Educ
> 
> ### include interaction
> #X<-as.matrix(dd2[,c(1,2,4)])
> ### don't include interaction
> X<-as.matrix(dd2[,c(1,2)])
> 
> ### create X and Y matricies
> X<-cbind(1,X)
> Y<-as.matrix(dd2[,3])
> 
> ### show X and Y matricies
> X
     Sex Educ
1  1   1    5
2  1   1    9
3  1   1   10
4  1   1   12
5  1   1   14
6  1   1   16
7  1   0    8
8  1   0   11
9  1   0   11
10 1   0   14
11 1   0   14
12 1   0   15
> Y
      [,1]
 [1,] 17.2
 [2,] 28.9
 [3,] 42.5
 [4,] 56.1
 [5,] 56.6
 [6,] 70.5
 [7,] 21.6
 [8,] 18.5
 [9,] 17.9
[10,] 31.8
[11,] 35.2
[12,] 39.5
> 
> ### find XtX
> XtX <- t(X)%*%X
> ### find XtY
> XtY <- t(X)%*%Y
> ### find (XtX)^-1
> XtX.inv <- solve(XtX)
> 
> ### show XtX and Xty.  These form the normal equations
> XtX
         Sex Educ
      12   6  139
Sex    6   6   66
Educ 139  66 1725
> XtY
       [,1]
      436.3
Sex   271.8
Educ 5468.4
> 
> ### find b, Y.hat and e
> b<-solve(XtX,XtY)
> Y.hat <- X%*%b
> e <- Y-Y.hat
> 
> ### show b, Y.hat and e
> b
           [,1]
     -24.965865
Sex   22.906316
Educ   4.305414
> Y.hat
        [,1]
1  19.467519
2  36.689173
3  40.994586
4  49.605414
5  58.216241
6  66.827068
7   9.477444
8  22.393684
9  22.393684
10 35.309925
11 35.309925
12 39.615338
> e
         [,1]
1  -2.2675188
2  -7.7891729
3   1.5054135
4   6.4945865
5  -1.6162406
6   3.6729323
7  12.1225564
8  -3.8936842
9  -4.4936842
10 -3.5099248
11 -0.1099248
12 -0.1153383
> 
> ### find MSE
> df<-length(Y)-ncol(X)
> MSE <- sum(e*e)/df
> 
> ### find t-statistic, p-values and CI's for b vector
> t.star<-b/sqrt(MSE*diag(XtX.inv))
> p.value<-2*pt(-abs(t.star),df)
> CI<-cbind(b-sqrt(MSE*diag(XtX.inv))*qt(.95,df),b+sqrt(MSE*diag(XtX.inv
+ ))*qt(.95,df))
> 
> ### show p-values and CI's
> p.value
             [,1]
     7.742132e-03
Sex  1.083812e-04
Educ 3.363518e-05
> CI
           [,1]       [,2]
     -38.384357 -11.547372
Sex   16.470146  29.342486
Educ   3.265501   5.345326
> 
> ### find hat matrix and I-H
> H<-X%*%XtX.inv%*%t(X)
> I.minus.H<-diag(length(Y))-H
> 
> ### show hat matrix and I-H and that their matrix product is 0
> H
              1           2           3           4            5           6           7
1   0.491478697  0.27493734  0.22080201  0.11253133  0.004260652 -0.10401003  0.22556391
2   0.274937343  0.20275689  0.18471178  0.14862155  0.112531328  0.07644110  0.07518797
3   0.220802005  0.18471178  0.17568922  0.15764411  0.139598997  0.12155388  0.03759398
4   0.112531328  0.14862155  0.15764411  0.17568922  0.193734336  0.21177945 -0.03759398
5   0.004260652  0.11253133  0.13959900  0.19373434  0.247869674  0.30200501 -0.11278195
6  -0.104010025  0.07644110  0.12155388  0.21177945  0.302005013  0.39223058 -0.18796992
7   0.225563910  0.07518797  0.03759398 -0.03759398 -0.112781955 -0.18796992  0.32330827
8   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632
9   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632
10 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436
11 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436
12 -0.153383459 -0.05112782 -0.02556391  0.02556391  0.076691729  0.12781955  0.06015038
             8           9          10          11          12
1   0.06315789  0.06315789 -0.09924812 -0.09924812 -0.15338346
2   0.02105263  0.02105263 -0.03308271 -0.03308271 -0.05112782
3   0.01052632  0.01052632 -0.01654135 -0.01654135 -0.02556391
4  -0.01052632 -0.01052632  0.01654135  0.01654135  0.02556391
5  -0.03157895 -0.03157895  0.04962406  0.04962406  0.07669173
6  -0.05263158 -0.05263158  0.08270677  0.08270677  0.12781955
7   0.21052632  0.21052632  0.09774436  0.09774436  0.06015038
8   0.17894737  0.17894737  0.14736842  0.14736842  0.13684211
9   0.17894737  0.17894737  0.14736842  0.14736842  0.13684211
10  0.14736842  0.14736842  0.19699248  0.19699248  0.21353383
11  0.14736842  0.14736842  0.19699248  0.19699248  0.21353383
12  0.13684211  0.13684211  0.21353383  0.21353383  0.23909774
> I.minus.H
              1           2           3           4            5           6           7
1   0.508521303 -0.27493734 -0.22080201 -0.11253133 -0.004260652  0.10401003 -0.22556391
2  -0.274937343  0.79724311 -0.18471178 -0.14862155 -0.112531328 -0.07644110 -0.07518797
3  -0.220802005 -0.18471178  0.82431078 -0.15764411 -0.139598997 -0.12155388 -0.03759398
4  -0.112531328 -0.14862155 -0.15764411  0.82431078 -0.193734336 -0.21177945  0.03759398
5  -0.004260652 -0.11253133 -0.13959900 -0.19373434  0.752130326 -0.30200501  0.11278195
6   0.104010025 -0.07644110 -0.12155388 -0.21177945 -0.302005013  0.60776942  0.18796992
7  -0.225563910 -0.07518797 -0.03759398  0.03759398  0.112781955  0.18796992  0.67669173
8  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632
9  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632
10  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436
11  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436
12  0.153383459  0.05112782  0.02556391 -0.02556391 -0.076691729 -0.12781955 -0.06015038
             8           9          10          11          12
1  -0.06315789 -0.06315789  0.09924812  0.09924812  0.15338346
2  -0.02105263 -0.02105263  0.03308271  0.03308271  0.05112782
3  -0.01052632 -0.01052632  0.01654135  0.01654135  0.02556391
4   0.01052632  0.01052632 -0.01654135 -0.01654135 -0.02556391
5   0.03157895  0.03157895 -0.04962406 -0.04962406 -0.07669173
6   0.05263158  0.05263158 -0.08270677 -0.08270677 -0.12781955
7  -0.21052632 -0.21052632 -0.09774436 -0.09774436 -0.06015038
8   0.82105263 -0.17894737 -0.14736842 -0.14736842 -0.13684211
9  -0.17894737  0.82105263 -0.14736842 -0.14736842 -0.13684211
10 -0.14736842 -0.14736842  0.80300752 -0.19699248 -0.21353383
11 -0.14736842 -0.14736842 -0.19699248  0.80300752 -0.21353383
12 -0.13684211 -0.13684211 -0.21353383 -0.21353383  0.76090226
> H%*%I.minus.H
               1            2            3            4            5             6
1   4.660968e-16 2.771670e-16 2.483865e-16 1.467862e-16 6.667335e-17 -1.998490e-17
2   2.581638e-16 2.043323e-16 2.113363e-16 2.018468e-16 1.694511e-16  1.364638e-16
3   2.099850e-16 2.080638e-16 2.275722e-16 1.884871e-16 1.839190e-16  1.972700e-16
4   1.260478e-16 1.624209e-16 1.974956e-16 2.156428e-16 2.474113e-16  2.712358e-16
5   4.222798e-17 1.272256e-16 1.787142e-16 2.351895e-16 3.208197e-16  3.524767e-16
6  -5.184180e-17 1.247540e-16 1.793306e-16 2.442290e-16 3.269276e-16  4.542570e-16
7   3.762258e-16 2.541965e-16 2.402612e-16 1.807484e-16 1.069078e-16  4.784466e-17
8   2.427681e-16 2.108265e-16 2.046576e-16 1.931053e-16 1.842373e-16  1.682665e-16
9   2.427681e-16 2.108265e-16 2.046576e-16 1.931053e-16 1.842373e-16  1.682665e-16
10  9.559275e-17 1.659194e-16 1.806454e-16 2.242287e-16 2.632223e-16  2.914454e-16
11  9.559275e-17 1.659194e-16 1.806454e-16 2.242287e-16 2.632223e-16  2.914454e-16
12  3.460299e-17 1.621992e-16 1.753748e-16 2.260875e-16 2.810608e-16  3.405225e-16
               7            8            9           10           11           12
1   4.086062e-16 3.115794e-16 3.115794e-16 1.929473e-16 1.929473e-16 1.220202e-16
2   2.659118e-16 2.129614e-16 2.129614e-16 1.737714e-16 1.737714e-16 1.678819e-16
3   2.004145e-16 1.937694e-16 1.955041e-16 1.835957e-16 1.835957e-16 1.767097e-16
4   1.075349e-16 1.478252e-16 1.460905e-16 1.825364e-16 1.825364e-16 1.923985e-16
5   1.105039e-17 9.886823e-17 9.886823e-17 1.818749e-16 1.818749e-16 2.240775e-16
6  -7.103599e-17 4.883314e-17 5.577204e-17 1.809245e-16 1.809245e-16 2.516840e-16
7   4.175159e-16 3.206816e-16 3.137427e-16 2.328892e-16 2.398281e-16 2.092985e-16
8   2.706643e-16 2.475742e-16 2.545131e-16 2.300186e-16 2.230797e-16 2.617196e-16
9   2.706643e-16 2.475742e-16 2.545131e-16 2.300186e-16 2.230797e-16 2.617196e-16
10  1.327894e-16 1.837994e-16 1.837994e-16 2.598155e-16 2.598155e-16 2.580266e-16
11  1.327894e-16 1.837994e-16 1.837994e-16 2.598155e-16 2.598155e-16 2.580266e-16
12  9.788313e-17 1.853410e-16 1.714632e-16 2.464764e-16 2.464764e-16 2.825431e-16
> 
> ### Cov(b,e) = E[be'] after a bit of manipulation
> ### we see that the numeric mean of each row is 0
> ### since each entry in a row is an iid random
> ### variable, taking the numeric mean of the row
> ### is a good estimate of the expected value of any
> ### entry in that row.
> mean(b[1]*e)
[1] -4.72955e-14
> mean(b[2]*e)
[1] 4.77762e-14
> mean(b[3]*e)
[1] 9.154547e-15
> ### We conclude numerically that E[be'] is the 0 matrix.
> ### I proved it symbolically, but am too lazy to type it
> ### into the Wiki.
> 

With Interaction Terms

> ### read input and make a copy ###
> dd <- read.table("data1.txt", header = T)
> dd2 <- dd
> 
> ### make required variables ###
> dd2$Sex <- ifelse(dd2$Sex == "M",1,0)
> dd2$sex.Educ <- dd2$Sex*dd2$Educ
> 
> ### include interaction
> X<-as.matrix(dd2[,c(1,2,4)])
> ### don't include interaction
> #X<-as.matrix(dd2[,c(1,2)])
> 
> ### create X and Y matricies
> X<-cbind(1,X)
> Y<-as.matrix(dd2[,3])
> 
> ### show X and Y matricies
> X
     Sex Educ sex.Educ
1  1   1    5        5
2  1   1    9        9
3  1   1   10       10
4  1   1   12       12
5  1   1   14       14
6  1   1   16       16
7  1   0    8        0
8  1   0   11        0
9  1   0   11        0
10 1   0   14        0
11 1   0   14        0
12 1   0   15        0
> Y
      [,1]
 [1,] 17.2
 [2,] 28.9
 [3,] 42.5
 [4,] 56.1
 [5,] 56.6
 [6,] 70.5
 [7,] 21.6
 [8,] 18.5
 [9,] 17.9
[10,] 31.8
[11,] 35.2
[12,] 39.5
> 
> ### find XtX
> XtX <- t(X)%*%X
> ### find XtY
> XtY <- t(X)%*%Y
> ### find (XtX)^-1
> XtX.inv <- solve(XtX)
> 
> ### show XtX and Xty.  These form the normal equations
> XtX
             Sex Educ sex.Educ
          12   6  139       66
Sex        6   6   66       66
Educ     139  66 1725      802
sex.Educ  66  66  802      802
> XtY
           [,1]
          436.3
Sex       271.8
Educ     5468.4
sex.Educ 3364.7
> 
> ### find b, Y.hat and e
> b<-solve(XtX,XtY)
> Y.hat <- X%*%b
> e <- Y-Y.hat
> 
> ### show b, Y.hat and e
> b
               [,1]
         -8.3090909
Sex      -0.6527512
Educ      2.9363636
sex.Educ  1.9965311
> Y.hat
       [,1]
1  15.70263
2  35.43421
3  40.36711
4  50.23289
5  60.09868
6  69.96447
7  15.18182
8  23.99091
9  23.99091
10 32.80000
11 32.80000
12 35.73636
> e
         [,1]
1   1.4973684
2  -6.5342105
3   2.1328947
4   5.8671053
5  -3.4986842
6   0.5355263
7   6.4181818
8  -5.4909091
9  -6.0909091
10 -1.0000000
11  2.4000000
12  3.7636364
> 
> ### find MSE
> df<-length(Y)-ncol(X)
> MSE <- sum(e*e)/df
> 
> ### find t-statistic, p-values and CI's for b vector
> t.star<-b/sqrt(MSE*diag(XtX.inv))
> p.value<-2*pt(-abs(t.star),df)
> CI<-cbind(b-sqrt(MSE*diag(XtX.inv))*qt(.95,df),b+sqrt(MSE*diag(XtX.inv
+ ))*qt(.95,df))
> 
> ### show p-values and CI's
> p.value
               [,1]
         0.47802864
Sex      0.96177830
Educ     0.01149265
sex.Educ 0.10357806
> CI
                [,1]      [,2]
         -29.0705049 12.452323
Sex      -25.2026159 23.897114
Educ       1.2624552  4.610272
sex.Educ  -0.0249054  4.017968
> 
> ### find hat matrix and I-H
> H<-X%*%XtX.inv%*%t(X)
> I.minus.H<-diag(length(Y))-H
> 
> ### show hat matrix and I-H and that their matrix product is 0
> H
               1            2            3            4             5             6
1   6.403509e-01 3.245614e-01 2.456140e-01 8.771930e-02 -7.017544e-02 -2.280702e-01
2   3.245614e-01 2.192982e-01 1.929825e-01 1.403509e-01  8.771930e-02  3.508772e-02
3   2.456140e-01 1.929825e-01 1.798246e-01 1.535088e-01  1.271930e-01  1.008772e-01
4   8.771930e-02 1.403509e-01 1.535088e-01 1.798246e-01  2.061404e-01  2.324561e-01
5  -7.017544e-02 8.771930e-02 1.271930e-01 2.061404e-01  2.850877e-01  3.640351e-01
6  -2.280702e-01 3.508772e-02 1.008772e-01 2.324561e-01  3.640351e-01  4.956140e-01
7   1.665335e-16 1.665335e-16 3.330669e-16 2.220446e-16  3.330669e-16  4.440892e-16
8   8.326673e-17 1.387779e-16 1.665335e-16 1.665335e-16  1.665335e-16  2.220446e-16
9   8.326673e-17 1.387779e-16 1.665335e-16 1.665335e-16  1.665335e-16  2.220446e-16
10  5.551115e-17 5.551115e-17 1.110223e-16 1.110223e-16  1.110223e-16  1.110223e-16
11  5.551115e-17 5.551115e-17 1.110223e-16 1.110223e-16  1.110223e-16  1.110223e-16
12  2.775558e-17 2.775558e-17 5.551115e-17 0.000000e+00 -5.551115e-17  0.000000e+00
               7             8             9            10            11            12
1   2.775558e-16  2.983724e-16  2.983724e-16  3.191891e-16  3.191891e-16  3.261280e-16
2   5.551115e-17  1.595946e-16  1.595946e-16  2.636780e-16  2.636780e-16  2.983724e-16
3  -3.330669e-16 -2.914335e-16 -2.914335e-16 -2.498002e-16 -2.498002e-16 -2.359224e-16
4  -6.661338e-16 -5.828671e-16 -5.828671e-16 -4.996004e-16 -4.996004e-16 -4.718448e-16
5  -1.110223e-16  1.387779e-17  1.387779e-17  1.387779e-16  1.387779e-16  1.804112e-16
6  -4.440892e-16 -2.775558e-16 -2.775558e-16 -1.110223e-16 -1.110223e-16 -5.551115e-17
7   6.650718e-01  3.062201e-01  3.062201e-01 -5.263158e-02 -5.263158e-02 -1.722488e-01
8   3.062201e-01  2.057416e-01  2.057416e-01  1.052632e-01  1.052632e-01  7.177033e-02
9   3.062201e-01  2.057416e-01  2.057416e-01  1.052632e-01  1.052632e-01  7.177033e-02
10 -5.263158e-02  1.052632e-01  1.052632e-01  2.631579e-01  2.631579e-01  3.157895e-01
11 -5.263158e-02  1.052632e-01  1.052632e-01  2.631579e-01  2.631579e-01  3.157895e-01
12 -1.722488e-01  7.177033e-02  7.177033e-02  3.157895e-01  3.157895e-01  3.971292e-01
> I.minus.H
               1             2             3             4             5             6
1   3.596491e-01 -3.245614e-01 -2.456140e-01 -8.771930e-02  7.017544e-02  2.280702e-01
2  -3.245614e-01  7.807018e-01 -1.929825e-01 -1.403509e-01 -8.771930e-02 -3.508772e-02
3  -2.456140e-01 -1.929825e-01  8.201754e-01 -1.535088e-01 -1.271930e-01 -1.008772e-01
4  -8.771930e-02 -1.403509e-01 -1.535088e-01  8.201754e-01 -2.061404e-01 -2.324561e-01
5   7.017544e-02 -8.771930e-02 -1.271930e-01 -2.061404e-01  7.149123e-01 -3.640351e-01
6   2.280702e-01 -3.508772e-02 -1.008772e-01 -2.324561e-01 -3.640351e-01  5.043860e-01
7  -1.665335e-16 -1.665335e-16 -3.330669e-16 -2.220446e-16 -3.330669e-16 -4.440892e-16
8  -8.326673e-17 -1.387779e-16 -1.665335e-16 -1.665335e-16 -1.665335e-16 -2.220446e-16
9  -8.326673e-17 -1.387779e-16 -1.665335e-16 -1.665335e-16 -1.665335e-16 -2.220446e-16
10 -5.551115e-17 -5.551115e-17 -1.110223e-16 -1.110223e-16 -1.110223e-16 -1.110223e-16
11 -5.551115e-17 -5.551115e-17 -1.110223e-16 -1.110223e-16 -1.110223e-16 -1.110223e-16
12 -2.775558e-17 -2.775558e-17 -5.551115e-17  0.000000e+00  5.551115e-17  0.000000e+00
               7             8             9            10            11            12
1  -2.775558e-16 -2.983724e-16 -2.983724e-16 -3.191891e-16 -3.191891e-16 -3.261280e-16
2  -5.551115e-17 -1.595946e-16 -1.595946e-16 -2.636780e-16 -2.636780e-16 -2.983724e-16
3   3.330669e-16  2.914335e-16  2.914335e-16  2.498002e-16  2.498002e-16  2.359224e-16
4   6.661338e-16  5.828671e-16  5.828671e-16  4.996004e-16  4.996004e-16  4.718448e-16
5   1.110223e-16 -1.387779e-17 -1.387779e-17 -1.387779e-16 -1.387779e-16 -1.804112e-16
6   4.440892e-16  2.775558e-16  2.775558e-16  1.110223e-16  1.110223e-16  5.551115e-17
7   3.349282e-01 -3.062201e-01 -3.062201e-01  5.263158e-02  5.263158e-02  1.722488e-01
8  -3.062201e-01  7.942584e-01 -2.057416e-01 -1.052632e-01 -1.052632e-01 -7.177033e-02
9  -3.062201e-01 -2.057416e-01  7.942584e-01 -1.052632e-01 -1.052632e-01 -7.177033e-02
10  5.263158e-02 -1.052632e-01 -1.052632e-01  7.368421e-01 -2.631579e-01 -3.157895e-01
11  5.263158e-02 -1.052632e-01 -1.052632e-01 -2.631579e-01  7.368421e-01 -3.157895e-01
12  1.722488e-01 -7.177033e-02 -7.177033e-02 -3.157895e-01 -3.157895e-01  6.028708e-01
> H%*%I.minus.H
               1             2             3             4             5             6
1   4.563509e-16  5.139169e-16  5.561466e-16  5.632464e-16  5.772834e-16  6.954073e-16
2   4.899044e-16  5.594520e-16  5.406846e-16  5.491120e-16  6.017754e-16  5.572935e-16
3   5.222449e-16  5.327109e-16  5.325466e-16  5.694382e-16  5.889827e-16  5.287315e-16
4   5.612982e-16  5.437867e-16  5.555604e-16  5.355281e-16  5.328447e-16  5.093446e-16
5   6.094233e-16  5.747500e-16  5.675595e-16  5.311371e-16  4.669659e-16  5.138170e-16
6   6.532860e-16  5.659366e-16  5.426093e-16  5.232224e-16  4.344398e-16  4.566931e-16
7  -1.219120e-16 -2.494018e-16 -2.233726e-16 -3.168651e-16 -3.359885e-16 -3.699858e-16
8  -1.072153e-16 -1.183707e-16 -1.776667e-16 -1.615977e-16 -2.135232e-16 -2.407475e-16
9  -1.072153e-16 -1.183707e-16 -1.776667e-16 -1.615977e-16 -2.135232e-16 -2.407475e-16
10 -4.869399e-17 -8.375367e-17 -6.914547e-17 -6.914547e-17 -5.745891e-17 -9.251859e-17
11 -4.869399e-17 -8.375367e-17 -6.914547e-17 -6.914547e-17 -5.745891e-17 -9.251859e-17
12 -4.590515e-17 -3.050015e-17 -1.828238e-17 -5.918533e-17 -6.184137e-17 -5.533408e-18
               7             8             9            10            11            12
1  -1.645857e-16 -1.824807e-16 -1.824807e-16 -2.003758e-16 -2.003758e-16 -2.063408e-16
2   8.083203e-17  1.472993e-17  1.472993e-17 -5.137216e-17 -5.137216e-17 -7.340619e-17
3   1.421865e-16  6.403260e-17  6.403260e-17 -1.412126e-17 -1.412126e-17 -4.017254e-17
4   2.648953e-16  1.626379e-16  1.626379e-16  6.038055e-17  6.038055e-17  2.629476e-17
5   3.876042e-16  2.612433e-16  2.612433e-16  1.348824e-16  1.348824e-16  9.276206e-17
6   5.103130e-16  3.598486e-16  3.598486e-16  2.093842e-16  2.093842e-16  1.592294e-16
7   2.660946e-15  1.279395e-15  1.275926e-15 -8.827777e-17 -8.827777e-17 -5.662991e-16
8   1.280731e-15  8.341009e-16  8.375703e-16  3.926743e-16  3.926743e-16  2.585111e-16
9   1.280731e-15  8.341009e-16  8.375703e-16  3.926743e-16  3.926743e-16  2.585111e-16
10 -9.770694e-17  3.843090e-16  3.843090e-16  8.732639e-16  8.732639e-16  1.066990e-15
11 -9.770694e-17  3.843090e-16  3.843090e-16  8.732639e-16  8.732639e-16  1.066990e-15
12 -5.552741e-16  2.252447e-16  2.252447e-16  1.033516e-15  1.061271e-15  1.286216e-15
> 
> ### Cov(b,e) = E[be'] after a bit of manipulation
> ### we see that the numeric mean of each row is 0
> ### since each entry in a row is an iid random
> ### variable, taking the numeric mean of the row
> ### is a good estimate of the expected value of any
> ### entry in that row.
> mean(b[1]*e)
[1] -2.723689e-14
> mean(b[2]*e)
[1] -2.169543e-15
> mean(b[3]*e)
[1] 1.002938e-14
> ### We conclude numerically that E[be'] is the 0 matrix.
> ### I proved it symbolically, but am too lazy to type it
> ### into the Wiki.
> 

Sbush 23:57, 21 Sep 2006 (EDT)

Male:1, Female:-1

Without Interaction Terms

> a<- read.csv("Q:\\R\\data.csv",header=T)
> a1<-a
###  X  Y  X'X  X'Y ###
> a1$Sex<-ifelse(a1$Sex=="M",1,-1)
> X<-as.matrix(a1[,c(1,2)])
> X<-cbind(1,X)
> X
     Sex Educ
1  1   1    5
2  1   1    9
3  1   1   10
4  1   1   12
5  1   1   14
6  1   1   16
7  1  -1    8
8  1  -1   11
9  1  -1   11
10 1  -1   14
11 1  -1   14
12 1  -1   15
> Y<-as.matrix(a1[,3])
> Y
     [,1]
[1,] 17.2
[2,] 28.9
[3,] 42.5
[4,] 56.1
[5,] 56.6
[6,] 70.5
[7,] 21.6
[8,] 18.5
[9,] 17.9
[10,] 31.8
[11,] 35.2
[12,] 39.5
> crossprod(X)
         Sex Educ
      12   0  139
Sex    0  12   -7
Educ 139  -7 1725
> crossprod(X,Y)
       [,1]
      436.3
Sex   107.3
Educ 5468.4
### b  Y.hat  e  H  I-H  H(I-H) ###
> b<-solve(crossprod(X),crossprod(X,Y))
> b
           [,1]
     -13.512707
Sex   11.453158
Educ   4.305414
> Y.hat<-X%*%b
> Y.hat
        [,1]
1  19.467519
2  36.689173
3  40.994586
4  49.605414
5  58.216241
6  66.827068
7  9.477444
8  22.393684
9  22.393684
10 35.309925
11 35.309925
12 39.615338
> e<-Y-Y.hat
> e
         [,1]
1  -2.2675188
2  -7.7891729
3   1.5054135
4   6.4945865
5  -1.6162406
6   3.6729323
7  12.1225564
8  -3.8936842
9  -4.4936842
10 -3.5099248
11 -0.1099248
12 -0.1153383
 
> H<-X%*%solve(crossprod(X))%*%t(X)
> H
              1           2           3           4            5           6           7           8           9
1   0.491478697  0.27493734  0.22080201  0.11253133  0.004260652 -0.10401003  0.22556391  0.06315789  0.06315789
2   0.274937343  0.20275689  0.18471178  0.14862155  0.112531328  0.07644110  0.07518797  0.02105263  0.02105263
3   0.220802005  0.18471178  0.17568922  0.15764411  0.139598997  0.12155388  0.03759398  0.01052632  0.01052632
4   0.112531328  0.14862155  0.15764411  0.17568922  0.193734336  0.21177945 -0.03759398 -0.01052632 -0.01052632
5   0.004260652  0.11253133  0.13959900  0.19373434  0.247869674  0.30200501 -0.11278195 -0.03157895 -0.03157895
6  -0.104010025  0.07644110  0.12155388  0.21177945  0.302005013  0.39223058 -0.18796992 -0.05263158 -0.05263158
7   0.225563910  0.07518797  0.03759398 -0.03759398 -0.112781955 -0.18796992  0.32330827  0.21052632  0.21052632
8   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632  0.17894737  0.17894737
9   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632  0.17894737  0.17894737
10 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436  0.14736842  0.14736842
11 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436  0.14736842  0.14736842
12 -0.153383459 -0.05112782 -0.02556391  0.02556391  0.076691729  0.12781955  0.06015038  0.13684211  0.13684211
            10          11          12
1  -0.09924812 -0.09924812 -0.15338346
2  -0.03308271 -0.03308271 -0.05112782
3  -0.01654135 -0.01654135 -0.02556391
4   0.01654135  0.01654135  0.02556391
5   0.04962406  0.04962406  0.07669173
6   0.08270677  0.08270677  0.12781955
7   0.09774436  0.09774436  0.06015038
8   0.14736842  0.14736842  0.13684211
9   0.14736842  0.14736842  0.13684211
10  0.19699248  0.19699248  0.21353383
11  0.19699248  0.19699248  0.21353383
12  0.21353383  0.21353383  0.23909774
> I<-diag(1,12,12)
> I-H
              1           2           3           4            5           6           7           8           9
1   0.508521303 -0.27493734 -0.22080201 -0.11253133 -0.004260652  0.10401003 -0.22556391 -0.06315789 -0.06315789
2  -0.274937343  0.79724311 -0.18471178 -0.14862155 -0.112531328 -0.07644110 -0.07518797 -0.02105263 -0.02105263
3  -0.220802005 -0.18471178  0.82431078 -0.15764411 -0.139598997 -0.12155388 -0.03759398 -0.01052632 -0.01052632
4  -0.112531328 -0.14862155 -0.15764411  0.82431078 -0.193734336 -0.21177945  0.03759398  0.01052632  0.01052632
5  -0.004260652 -0.11253133 -0.13959900 -0.19373434  0.752130326 -0.30200501  0.11278195  0.03157895  0.03157895
6   0.104010025 -0.07644110 -0.12155388 -0.21177945 -0.302005013  0.60776942  0.18796992  0.05263158  0.05263158
7  -0.225563910 -0.07518797 -0.03759398  0.03759398  0.112781955  0.18796992  0.67669173 -0.21052632 -0.21052632
8  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632  0.82105263 -0.17894737
9  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632 -0.17894737  0.82105263
10  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436 -0.14736842 -0.14736842
11  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436 -0.14736842 -0.14736842
12  0.153383459  0.05112782  0.02556391 -0.02556391 -0.076691729 -0.12781955 -0.06015038 -0.13684211 -0.13684211
            10          11          12
1   0.09924812  0.09924812  0.15338346
2   0.03308271  0.03308271  0.05112782
3   0.01654135  0.01654135  0.02556391
4  -0.01654135 -0.01654135 -0.02556391
5  -0.04962406 -0.04962406 -0.07669173
6  -0.08270677 -0.08270677 -0.12781955
7  -0.09774436 -0.09774436 -0.06015038
8  -0.14736842 -0.14736842 -0.13684211
9  -0.14736842 -0.14736842 -0.13684211
10  0.80300752 -0.19699248 -0.21353383
11 -0.19699248  0.80300752 -0.21353383
12 -0.21353383 -0.21353383  0.76090226
> H%*%(I-H)
              1             2             3             4             5             6             7             8
1   1.461030e-16  2.400745e-17 -8.527928e-17 -1.798192e-16 -3.047167e-16 -3.983901e-16  1.004852e-16 -5.717472e-17
2  -3.200938e-18 -1.141525e-16 -1.666990e-16 -1.893928e-16 -2.537198e-16 -3.093733e-16 -3.402976e-17 -1.258289e-16
3  -6.630023e-17 -1.332868e-16 -1.583176e-16 -2.161857e-16 -2.519360e-16 -2.876863e-16 -7.791930e-17 -1.469799e-16
4  -1.558943e-16 -1.662922e-16 -1.940008e-16 -2.158077e-16 -2.415179e-16 -2.297145e-16 -1.538668e-16 -1.773221e-16
5  -2.517365e-16 -2.210404e-16 -2.271111e-16 -2.106300e-16 -1.928478e-16 -2.180000e-16 -2.160586e-16 -2.109510e-16
6  -3.188368e-16 -2.714126e-16 -2.749473e-16 -2.260716e-16 -2.118904e-16 -1.478360e-16 -3.011541e-16 -2.441234e-16
7   2.441412e-16  1.276777e-16  1.080988e-16  2.644013e-17 -2.052403e-17 -8.917224e-17  2.066623e-16  1.258242e-16
8   9.980589e-17  5.649583e-17  4.590559e-17  2.299038e-17  1.809262e-18 -2.284109e-17  1.014729e-16  7.048669e-17
9   9.980589e-17  5.649583e-17  4.590559e-17  2.299038e-17  1.809262e-18 -2.284109e-17  1.014729e-16  7.048669e-17
10 -4.468607e-17 -1.488322e-17 -9.906050e-19  1.638670e-17  3.723387e-17  6.154880e-17  3.786237e-18  4.051359e-17
11 -4.468607e-17 -1.488322e-17 -9.906050e-19  1.638670e-17  3.723387e-17  6.154880e-17  3.786237e-18  4.051359e-17
12 -8.643802e-17 -3.904483e-17 -2.287836e-17  1.465833e-17  5.046114e-17  7.932294e-17 -4.443196e-17  1.455203e-17
               9            10            11            12
1  -5.023583e-17 -2.381721e-16 -2.381721e-16 -2.765325e-16
2  -1.258289e-16 -2.134091e-16 -2.168786e-16 -2.342182e-16
3  -1.469799e-16 -2.104619e-16 -2.121966e-16 -2.217922e-16 
4  -1.773221e-16 -1.992005e-16 -1.992005e-16 -1.989765e-16
5  -2.109510e-16 -1.993780e-16 -1.993780e-16 -1.914735e-16
6  -2.510623e-16 -1.835453e-16 -1.766064e-16 -1.887935e-16
7   1.119464e-16  2.280297e-17  1.933353e-17  1.101482e-17
8   7.742559e-17  4.329863e-17  5.370697e-17  2.498408e-17
9   7.742559e-17  4.329863e-17  5.370697e-17  2.498408e-17
10  2.663580e-17  5.888234e-17  5.194345e-17  7.425430e-17
11  2.663580e-17  5.888234e-17  5.194345e-17  7.425430e-17
12  1.455203e-17  4.721023e-17  4.721023e-17  9.973305e-17
### the above matrix is very close to 0 matrix in numerical, indeed H and I-H are orthogonal###
           
>  ### find t-statistic, p-values and CI's for b vector ###
> df<-length(Y)-ncol(X)
> df<-length(Y)-ncol(X)
> MSE <- sum(e*e)/df
> t.star<-b/sqrt(MSE*diag(solve(crossprod(X))))
> t<-b/sqrt(MSE*diag(solve(crossprod(X))))
> t
          [,1]
     -1.989050
Sex   6.524045
Educ  7.589401
  > CI<-cbind(b-sqrt(MSE*diag(solve(crossprod(X))))*qt(.95,df),b+sqrt(MSE*diag((solve(crossprod(X))*qt(.95,df))
+ )
+ )
+ )

> p<-2*pt(-abs(t.star),df)
> p
             [,1]
     7.791577e-02
Sex  1.083812e-04
Educ 3.363518e-05
### the p-value of b_0 is 7.791577e-02>0.05, then we accept the null hopothesis b_0=0 ###
### the p-value of b_1 is 1.083812e-04<0.05, then we refuse the null hopothesis b_1=0 ###
### the p-value of b_2 is 3.363518e-05<0.05, then we refuse the null hopothesis b_2=0 ###
> CI
           [,1]      [,2]
     -25.966049 -4.314751
Sex    8.235073 13.830014
Educ   3.265501  5.073486
> Cov<-mean(b%*%t(e))-mean(b)%*%t(mean(e))
> Cov
             [,1]
[1,] -2.722945e-16
### the Cov(b,e) is very close to 0, then numercalily we can conclude that b and  e are independent ###

With Interaction Terms

> a<- read.csv("Q:\\R\\data.csv",header=T)
> a1<-a
>###  X  Y  X'X  X'Y ### 
> a1$Sex<-ifelse(a1$Sex=="M",1,-1)
> a1$Sex.Educ<-a1$Sex*a1$Educ
> X<-as.matrix(a1[,c(1,2,4)])
> X<-cbind(1,X)
> X
     Sex Educ Sex.Educ
1  1   1    5        5
2  1   1    9        9
3  1   1   10       10
4  1   1   12       12
5  1   1   14       14
6  1   1   16       16
7  1  -1    8       -8
8  1  -1   11      -11
9  1  -1   11      -11
10 1  -1   14      -14
11 1  -1   14      -14
12 1  -1   15      -15
> Y<-as.matrix(a1[,3])
> Y
     [,1]
[1,] 17.2
[2,] 28.9
[3,] 42.5
[4,] 56.1
[5,] 56.6
[6,] 70.5
[7,] 21.6
[8,] 18.5
[9,] 17.9
[10,] 31.8
[11,] 35.2
[12,] 39.5
> crossprod(X)
             Sex Educ Sex.Educ
          12   0  139       -7
Sex        0   12  -7      139
Educ     139  -7  1725     -121
Sex.Educ  -7  139 -121     1725
>  crossprod(X,Y)
           [,1]
          436.3
Sex       107.3
Educ     5468.4
Sex.Educ 1261.0
### b  Y.hat  e  H  I-H  H(I-H) ###
> b<-solve(crossprod(X),crossprod(X,Y))
> b
               [,1]
         -8.6354665
Sex      -0.3263756
Educ      3.9346292
Sex.Educ  0.9982656
> Y.hat<-X%*%b
> Y.hat
       [,1]
1  15.70263
2  35.43421
3  40.36711
4  50.23289
5  60.09868
6  69.96447
7  15.18182
8  23.99091
9  23.99091
10 32.80000
11 32.80000
12 35.73636
> e<-Y-Y.hat
> e
         [,1]
1   1.4973684
2  -6.5342105
3   2.1328947
4   5.8671053
5  -3.4986842
6   0.5355263
7   6.4181818
8  -5.4909091
9  -6.0909091
10 -1.0000000
11  2.4000000
12  3.7636364 
> H<-X%*%solve(crossprod(X))%*%t(X)
> H
               1            2            3            4             5             6             7             8
1   6.403509e-01 3.245614e-01 2.456140e-01 8.771930e-02 -7.017544e-02 -2.280702e-01 -5.551115e-16 -1.526557e-16
2   3.245614e-01 2.192982e-01 1.929825e-01 1.403509e-01  8.771930e-02  3.508772e-02 -2.081668e-16 -5.724587e-17
3   2.456140e-01 1.929825e-01 1.798246e-01 1.535088e-01  1.271930e-01  1.008772e-01  6.938894e-18  9.367507e-17
4   8.771930e-02 1.403509e-01 1.535088e-01 1.798246e-01  2.061404e-01  2.324561e-01  9.714451e-17  4.857226e-17
5  -7.017544e-02 8.771930e-02 1.271930e-01 2.061404e-01  2.850877e-01  3.640351e-01  3.053113e-16  1.144917e-16
6  -2.280702e-01 3.508772e-02 1.008772e-01 2.324561e-01  3.640351e-01  4.956140e-01  4.996004e-16  1.942890e-16
7   6.938894e-17 4.579670e-16 5.828671e-16 7.216450e-16  8.604228e-16  1.110223e-15  6.650718e-01  3.062201e-01
8   2.081668e-16 2.914335e-16 3.053113e-16 3.608225e-16  4.163336e-16  4.440892e-16  3.062201e-01  2.057416e-01
9   2.081668e-16 2.914335e-16 3.053113e-16 3.608225e-16  4.163336e-16  4.440892e-16  3.062201e-01  2.057416e-01
10  5.238865e-16 3.989864e-16 3.816392e-16 2.914335e-16  2.012279e-16  1.665335e-16 -5.263158e-02  1.052632e-01
11  5.238865e-16 3.989864e-16 3.816392e-16 2.914335e-16  2.012279e-16  1.665335e-16 -5.263158e-02  1.052632e-01
12  7.979728e-16 6.036838e-16 5.967449e-16 4.718448e-16  4.579670e-16  3.330669e-16 -1.722488e-01  7.177033e-02
               9            10            11            12
1  -1.526557e-16  2.498002e-16  2.498002e-16  4.024558e-16
2  -5.724587e-17  9.367507e-17  9.367507e-17  1.439820e-16
3   9.367507e-17  1.734723e-16  1.734723e-16  2.046974e-16
4    4.857226e-17  0.000000e+00  0.000000e+00 -6.938894e-18
5   1.144917e-16 -7.632783e-17 -7.632783e-17 -1.214306e-16
6   1.942890e-16 -1.110223e-16 -1.110223e-16 -2.498002e-16
7   3.062201e-01 -5.263158e-02 -5.263158e-02 -1.722488e-01
8   2.057416e-01  1.052632e-01  1.052632e-01  7.177033e-02
9   2.057416e-01  1.052632e-01  1.052632e-01  7.177033e-02
10  1.052632e-01  2.631579e-01  2.631579e-01  3.157895e-01
11  1.052632e-01  2.631579e-01  2.631579e-01  3.157895e-01
12  7.177033e-02  3.157895e-01  3.157895e-01  3.971292e-01 
> I<-diag(1,12,12)
> I-H
               1             2             3             4             5             6             7             8
1   3.596491e-01 -3.245614e-01 -2.456140e-01 -8.771930e-02  7.017544e-02  2.280702e-01  5.551115e-16  1.526557e-16
2  -3.245614e-01  7.807018e-01 -1.929825e-01 -1.403509e-01 -8.771930e-02 -3.508772e-02  2.081668e-16  5.724587e-17
3  -2.456140e-01 -1.929825e-01  8.201754e-01 -1.535088e-01 -1.271930e-01 -1.008772e-01 -6.938894e-18 -9.367507e-17
4  -8.771930e-02 -1.403509e-01 -1.535088e-01  8.201754e-01 -2.061404e-01 -2.324561e-01 -9.714451e-17 -4.857226e-17
5   7.017544e-02 -8.771930e-02 -1.271930e-01 -2.061404e-01  7.149123e-01 -3.640351e-01 -3.053113e-16 -1.144917e-16
6   2.280702e-01 -3.508772e-02 -1.008772e-01 -2.324561e-01 -3.640351e-01  5.043860e-01 -4.996004e-16 -1.942890e-16
7  -6.938894e-17 -4.579670e-16 -5.828671e-16 -7.216450e-16 -8.604228e-16 -1.110223e-15  3.349282e-01 -3.062201e-01
8  -2.081668e-16 -2.914335e-16 -3.053113e-16 -3.608225e-16 -4.163336e-16 -4.440892e-16 -3.062201e-01  7.942584e-01
9  -2.081668e-16 -2.914335e-16 -3.053113e-16 -3.608225e-16 -4.163336e-16 -4.440892e-16 -3.062201e-01 -2.057416e-01
10 -5.238865e-16 -3.989864e-16 -3.816392e-16 -2.914335e-16 -2.012279e-16 -1.665335e-16  5.263158e-02 -1.052632e-01
11 -5.238865e-16 -3.989864e-16 -3.816392e-16 -2.914335e-16 -2.012279e-16 -1.665335e-16  5.263158e-02 -1.052632e-01
12 -7.979728e-16 -6.036838e-16 -5.967449e-16 -4.718448e-16 -4.579670e-16 -3.330669e-16  1.722488e-01 -7.177033e-02
               9            10            11            12
1   1.526557e-16 -2.498002e-16 -2.498002e-16 -4.024558e-16
2   5.724587e-17 -9.367507e-17 -9.367507e-17 -1.439820e-16
3  -9.367507e-17 -1.734723e-16 -1.734723e-16 -2.046974e-16
4  -4.857226e-17  0.000000e+00  0.000000e+00  6.938894e-18
5  -1.144917e-16  7.632783e-17  7.632783e-17  1.214306e-16
6  -1.942890e-16  1.110223e-16  1.110223e-16  2.498002e-16
7  -3.062201e-01  5.263158e-02  5.263158e-02  1.722488e-01
8  -2.057416e-01 -1.052632e-01 -1.052632e-01 -7.177033e-02
9   7.942584e-01 -1.052632e-01 -1.052632e-01 -7.177033e-02
10 -1.052632e-01  7.368421e-01 -2.631579e-01 -3.157895e-01
11 -1.052632e-01 -2.631579e-01  7.368421e-01 -3.157895e-01
12 -7.177033e-02 -3.157895e-01 -3.157895e-01  6.028708e-01
> H%*%(I-H)
              1             2             3             4             5             6             7             8
1  -2.387549e-16 -3.618144e-16 -3.855948e-16 -4.766153e-16 -5.051840e-16 -5.747762e-16  5.513599e-16  1.400824e-16
2  -1.031280e-16 -2.563854e-16 -3.174709e-16 -3.649348e-16 -4.245414e-16 -4.901136e-16  1.665335e-16  2.034496e-17
3  -8.576039e-17 -2.251274e-16 -2.814648e-16 -3.236089e-16 -4.021814e-16 -4.953042e-16  6.959644e-17 -8.128577e-18
4  -5.498938e-18 -1.987055e-16 -1.939604e-16 -2.975017e-16 -3.698214e-16 -4.495644e-16 -1.181023e-16 -7.111813e-17
5   5.151316e-17 -1.306430e-16 -1.679632e-16 -2.590430e-16 -3.362382e-16 -4.033368e-16 -3.073283e-16 -1.323149e-16
6   1.160706e-16 -4.078972e-17 -1.218203e-16 -2.141299e-16 -2.925584e-16 -3.711089e-16 -5.077096e-16 -1.888635e-16
7   2.401499e-18 -3.302670e-16 -3.912562e-16 -5.923846e-16 -7.743894e-16 -9.154913e-16 -7.647098e-16 -7.656204e-16
8  -2.703291e-16 -3.891867e-16 -4.373938e-16 -4.638213e-16 -4.982170e-16 -5.866630e-16 -5.660103e-16 -6.778601e-16
9  -2.703291e-16 -3.891867e-16 -4.373938e-16 -4.638213e-16 -4.982170e-16 -5.866630e-16 -5.660103e-16 -6.778601e-16
10 -5.761717e-16 -4.344721e-16 -4.026993e-16 -3.420753e-16 -3.165110e-16 -2.179056e-16 -3.990000e-16 -5.823894e-16
11 -5.761717e-16 -4.344721e-16 -4.026993e-16 -3.420753e-16 -3.165110e-16 -2.179056e-16 -3.990000e-16 -5.823894e-16
12 -6.545026e-16 -4.756186e-16 -4.087861e-16 -3.301010e-16 -1.844838e-16 -1.177509e-16 -3.119182e-16 -5.227447e-16
               9            10            11            12
1   1.400824e-16 -2.694908e-16 -2.694908e-16 -4.084499e-16
2   2.034496e-17 -1.245045e-16 -1.245045e-16 -1.813090e-16
3  -8.128577e-18 -9.154471e-17 -9.154471e-17 -1.201414e-16
4  -7.111813e-17 -2.306878e-17 -2.306878e-17  2.990807e-18
5  -1.323149e-16  4.358112e-17  4.358112e-17  1.277166e-16
6  -1.888635e-16  1.306825e-16  1.306825e-16  2.133987e-16
7  -7.794982e-16 -7.873476e-16 -7.942865e-16 -8.074663e-16
8  -6.639824e-16 -7.619536e-16 -7.619536e-16 -7.986911e-16
9  -6.639824e-16 -7.619536e-16 -7.619536e-16 -7.986911e-16
10 -5.685116e-16 -7.796566e-16 -7.657788e-16 -7.984878e-16
11 -5.685116e-16 -7.796566e-16 -7.657788e-16 -7.984878e-16
12 -5.227447e-16 -7.543879e-16 -7.543879e-16 -8.333042e-16
              
>  ### find t-statistic, p-values and CI's for b vector ###
> df<-length(Y)-ncol(X)
> MSE <- crossprod(e)/df
> t.star<-b/sqrt(MSE*diag(solve(crossprod(X))))
> t<-b/sqrt(MSE*diag(solve(crossprod(X))))
> t
                [,1]
         -1.30819986
Sex      -0.04944313
Educ      7.23904211
Sex.Educ  1.83663721
 > CI<-cbind(b-sqrt(MSE*diag(solve(crossprod(X))))*qt(.95,df),b+sqrt(MSE*diag((solve(crossprod(X))*qt(.95,df))
+ )))
> p<-2*pt(-abs(t.star),df)
> p
                 [,1]
         2.271398e-01
Sex      9.617783e-01
Educ     8.899721e-05
Sex.Educ 1.035781e-01
### the p-value of b_0 is 2.271398e-01>0.05, then we fail to reject the null hopothesis b_0=0 ###
### the p-value of b_1 is 9.617783e-01>0.05,  then we accept the null hopothesis b_1=0 ###
### the p-value of b_2 is 8.899721e-05<0.05,  then we fail to reject the null hopothesis b_2=0 ###
### the p-value of b_3 is 1.035781e-01>0.05,  then we fail to reject the null hopothesis b_3=0 ###
>CI
                [,1]      [,2]
         -20.9103989 0.3660443
Sex      -12.6013080 8.6751352
Educ       2.9239109 4.6758138
Sex.Educ  -0.0124527 1.7394502
> Cov<-mean(b%*%t(e))-mean(b)%*%t(mean(e))
> Cov
             [,1]
[1,] -1.128727e-16
### the Cov(b,e) is close to 0, then numercalily we can conclude that b and  e are independent ###

fang chang Sep.22

Male:2, Female:1

Without Interaction Terms

Y vector and X matrix

>  dd<-read.csv("Q:\\6630a1.csv",header=T)
>  dd2<-dd
>  y<-as.matrix(dd2[,4])
>  dd2$x1<-ifelse(dd2$Sex=="M",2,1)
>  dd2$x2<-dd2$Educ
>  xx<-as.matrix(dd2[,c(5,6)])
>  x<-cbind(1,xx)
>  dd2
    X Sex Educ  Sal x1 x2
1   1   M    5 17.2  2  5
2   2   M    9 28.9  2  9
3   3   M   10 42.5  2 10
4   4   M   12 56.1  2 12
5   5   M   14 56.6  2 14
6   6   M   16 70.5  2 16
7   7   F    8 21.6  1  8
8   8   F   11 18.5  1 11
9   9   F   11 17.9  1 11
10 10   F   14 31.8  1 14
11 11   F   14 35.2  1 14
12 12   F   15 39.5  1 15
> y
##here is y vector##
      [,1]
 [1,] 17.2
 [2,] 28.9
 [3,] 42.5
 [4,] 56.1
 [5,] 56.6
 [6,] 70.5
 [7,] 21.6
 [8,] 18.5
 [9,] 17.9
[10,] 31.8
[11,] 35.2
12,] 39.5
> x
##here is x matrix##
     x1 x2
1  1  2  5
2  1  2  9
3  1  2 10
4  1  2 12
5  1  2 14
6  1  2 16
7  1  1  8
8  1  1 11
9  1  1 11
10 1  1 14
11 1  1 14
12 1  1 15

X'X and X'Y

> xtx<-t(x)%*%x
> xtx
##here is x'x ##
        x1   x2
    12  18  139
x1  18  30  205
x2 139 205 1725
> xty<- t(x)%*%y
> xty
##here is x'y ##
     [,1]
    436.3
x1  708.1
x2 5468.4

\hat\beta

> xtxinv<-solve(xtx)
> xtxinv
                       x1           x2
    2.4373434 -0.64035088 -0.120300752
x1 -0.6403509  0.34561404  0.010526316
x2 -0.1203008  0.01052632  0.009022556
> b<-xtxinv%*%xty
> b
OR
> b<-solve(xtx,xty)
> b
##here is b hat ##
         [,1]
   -47.872180
x1  22.906316
x2   4.305414

Normal Equation

> nn<-xtx%*%b
> nn
     [,1]
    436.3
x1  708.1
x2 5468.4
##so we verify the normal equation x'xb=x'y here##

residual e and \hat y

> yhat<-x%*%b
> yhat
##here is yhat ##
        [,1]
1  19.467519
2  36.689173
3  40.994586
4  49.605414
5  58.216241
6  66.827068
7   9.477444
8  22.393684
9  22.393684
10 35.309925
11 35.309925
12 39.615338
> e<-y-yhat
> e
##here is residual e ##
         [,1]
1  -2.2675188
2  -7.7891729
3   1.5054135
4   6.4945865
5  -1.6162406
6   3.6729323
7  12.1225564
8  -3.8936842
9  -4.4936842
10 -3.5099248
11 -0.1099248
12 -0.1153383

H and (I-H)

> n<-length(y)
> n
[1] 12
> k<-ncol(x)
> k
[1] 3
> i<-diag(n)
> i
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,]    1    0    0    0    0    0    0    0    0     0     0     0
 [2,]    0    1    0    0    0    0    0    0    0     0     0     0
 [3,]    0    0    1    0    0    0    0    0    0     0     0     0
 [4,]    0    0    0    1    0    0    0    0    0     0     0     0
 [5,]    0    0    0    0    1    0    0    0    0     0     0     0
 [6,]    0    0    0    0    0    1    0    0    0     0     0     0
 [7,]    0    0    0    0    0    0    1    0    0     0     0     0
 [8,]    0    0    0    0    0    0    0    1    0     0     0     0
 [9,]    0    0    0    0    0    0    0    0    1     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     1     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     1     0
[12,]    0    0    0    0    0    0    0    0    0     0     0     1
> h<-x%*%xtxinv%*%t(x)
> h
##here is H matrix##
              1           2           3           4            5           6           7           8           9           10          11          12
1   0.491478697  0.27493734  0.22080201  0.11253133  0.004260652 -0.10401003  0.22556391  0.06315789  0.06315789 -0.09924812 -0.09924812 -0.15338346
2   0.274937343  0.20275689  0.18471178  0.14862155  0.112531328  0.07644110  0.07518797  0.02105263  0.02105263 -0.03308271 -0.03308271 -0.05112782
3   0.220802005  0.18471178  0.17568922  0.15764411  0.139598997  0.12155388  0.03759398  0.01052632  0.01052632 -0.01654135 -0.01654135 -0.02556391
4   0.112531328  0.14862155  0.15764411  0.17568922  0.193734336  0.21177945 -0.03759398 -0.01052632 -0.01052632  0.01654135  0.01654135  0.02556391
5   0.004260652  0.11253133  0.13959900  0.19373434  0.247869674  0.30200501 -0.11278195 -0.03157895 -0.03157895  0.04962406  0.04962406  0.07669173
6  -0.104010025  0.07644110  0.12155388  0.21177945  0.302005013  0.39223058 -0.18796992 -0.05263158 -0.05263158  0.08270677  0.08270677  0.12781955
7   0.225563910  0.07518797  0.03759398 -0.03759398 -0.112781955 -0.18796992  0.32330827  0.21052632  0.21052632  0.09774436  0.09774436  0.06015038
8   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632  0.17894737  0.17894737  0.14736842  0.14736842  0.13684211
9   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632  0.17894737  0.17894737  0.14736842  0.14736842  0.13684211
10 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436  0.14736842  0.14736842  0.19699248  0.19699248  0.21353383
11 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436  0.14736842  0.14736842  0.19699248  0.19699248  0.21353383
12 -0.153383459 -0.05112782 -0.02556391  0.02556391  0.076691729  0.12781955  0.06015038  0.13684211  0.13684211  0.21353383  0.21353383  0.23909774
> he<-i-h
> he
##here is I-H matrix##
              1           2           3           4            5           6           7           8           9           10          11          12
1   0.508521303 -0.27493734 -0.22080201 -0.11253133 -0.004260652  0.10401003 -0.22556391 -0.06315789 -0.06315789  0.09924812  0.09924812  0.15338346
2  -0.274937343  0.79724311 -0.18471178 -0.14862155 -0.112531328 -0.07644110 -0.07518797 -0.02105263 -0.02105263  0.03308271  0.03308271  0.05112782
3  -0.220802005 -0.18471178  0.82431078 -0.15764411 -0.139598997 -0.12155388 -0.03759398 -0.01052632 -0.01052632  0.01654135  0.01654135  0.02556391
4  -0.112531328 -0.14862155 -0.15764411  0.82431078 -0.193734336 -0.21177945  0.03759398  0.01052632  0.01052632 -0.01654135 -0.01654135 -0.02556391
5  -0.004260652 -0.11253133 -0.13959900 -0.19373434  0.752130326 -0.30200501  0.11278195  0.03157895  0.03157895 -0.04962406 -0.04962406 -0.07669173
6   0.104010025 -0.07644110 -0.12155388 -0.21177945 -0.302005013  0.60776942  0.18796992  0.05263158  0.05263158 -0.08270677 -0.08270677 -0.12781955
7  -0.225563910 -0.07518797 -0.03759398  0.03759398  0.112781955  0.18796992  0.67669173 -0.21052632 -0.21052632 -0.09774436 -0.09774436 -0.06015038
8  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632  0.82105263 -0.17894737 -0.14736842 -0.14736842 -0.13684211
9  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632 -0.17894737  0.82105263 -0.14736842 -0.14736842 -0.13684211
10  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436 -0.14736842 -0.14736842  0.80300752 -0.19699248 -0.21353383
11  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436 -0.14736842 -0.14736842 -0.19699248  0.80300752 -0.21353383
12  0.153383459  0.05112782  0.02556391 -0.02556391 -0.076691729 -0.12781955 -0.06015038 -0.13684211 -0.13684211 -0.21353383 -0.21353383  0.76090226
> h%*%h
##here is HH matrix##
              1           2           3           4            5           6           7           8           9          10          11          12
1   0.491478697  0.27493734  0.22080201  0.11253133  0.004260652 -0.10401003  0.22556391  0.06315789  0.06315789 -0.09924812 -0.09924812 -0.15338346
2   0.274937343  0.20275689  0.18471178  0.14862155  0.112531328  0.07644110  0.07518797  0.02105263  0.02105263 -0.03308271 -0.03308271 -0.05112782
3   0.220802005  0.18471178  0.17568922  0.15764411  0.139598997  0.12155388  0.03759398  0.01052632  0.01052632 -0.01654135 -0.01654135 -0.02556391
4   0.112531328  0.14862155  0.15764411  0.17568922  0.193734336  0.21177945 -0.03759398 -0.01052632 -0.01052632  0.01654135  0.01654135  0.02556391
5   0.004260652  0.11253133  0.13959900  0.19373434  0.247869674  0.30200501 -0.11278195 -0.03157895 -0.03157895  0.04962406  0.04962406  0.07669173
6  -0.104010025  0.07644110  0.12155388  0.21177945  0.302005013  0.39223058 -0.18796992 -0.05263158 -0.05263158  0.08270677  0.08270677  0.12781955
7   0.225563910  0.07518797  0.03759398 -0.03759398 -0.112781955 -0.18796992  0.32330827  0.21052632  0.21052632  0.09774436  0.09774436  0.06015038
8   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632  0.17894737  0.17894737  0.14736842  0.14736842  0.13684211
9   0.063157895  0.02105263  0.01052632 -0.01052632 -0.031578947 -0.05263158  0.21052632  0.17894737  0.17894737  0.14736842  0.14736842  0.13684211
10 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436  0.14736842  0.14736842  0.19699248  0.19699248  0.21353383
11 -0.099248120 -0.03308271 -0.01654135  0.01654135  0.049624060  0.08270677  0.09774436  0.14736842  0.14736842  0.19699248  0.19699248  0.21353383
12 -0.153383459 -0.05112782 -0.02556391  0.02556391  0.076691729  0.12781955  0.06015038  0.13684211  0.13684211  0.21353383  0.21353383  0.23909774
##By HH=H, we verify that H is idenpotent matrix##
> he%*%he
##here is (I-H)(I-H) matrix##
              1           2           3           4            5           6           7           8           9           10          11          12
1   0.508521303 -0.27493734 -0.22080201 -0.11253133 -0.004260652  0.10401003 -0.22556391 -0.06315789 -0.06315789  0.09924812  0.09924812  0.15338346
2  -0.274937343  0.79724311 -0.18471178 -0.14862155 -0.112531328 -0.07644110 -0.07518797 -0.02105263 -0.02105263  0.03308271  0.03308271  0.05112782
3  -0.220802005 -0.18471178  0.82431078 -0.15764411 -0.139598997 -0.12155388 -0.03759398 -0.01052632 -0.01052632  0.01654135  0.01654135  0.02556391
4  -0.112531328 -0.14862155 -0.15764411  0.82431078 -0.193734336 -0.21177945  0.03759398  0.01052632  0.01052632 -0.01654135 -0.01654135 -0.02556391
5  -0.004260652 -0.11253133 -0.13959900 -0.19373434  0.752130326 -0.30200501  0.11278195  0.03157895  0.03157895 -0.04962406 -0.04962406 -0.07669173
6   0.104010025 -0.07644110 -0.12155388 -0.21177945 -0.302005013  0.60776942  0.18796992  0.05263158  0.05263158 -0.08270677 -0.08270677 -0.12781955
7  -0.225563910 -0.07518797 -0.03759398  0.03759398  0.112781955  0.18796992  0.67669173 -0.21052632 -0.21052632 -0.09774436 -0.09774436 -0.06015038
8  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632  0.82105263 -0.17894737 -0.14736842 -0.14736842 -0.13684211
9  -0.063157895 -0.02105263 -0.01052632  0.01052632  0.031578947  0.05263158 -0.21052632 -0.17894737  0.82105263 -0.14736842 -0.14736842 -0.13684211
10  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436 -0.14736842 -0.14736842  0.80300752 -0.19699248 -0.21353383
11  0.099248120  0.03308271  0.01654135 -0.01654135 -0.049624060 -0.08270677 -0.09774436 -0.14736842 -0.14736842 -0.19699248  0.80300752 -0.21353383
12  0.153383459  0.05112782  0.02556391 -0.02556391 -0.076691729 -0.12781955 -0.06015038 -0.13684211 -0.13684211 -0.21353383 -0.21353383  0.76090226
##By (I-H)(I-H)=I-H, we verify that I-H is idenpotent matrix##
> h%*%he
##here is H(I-H) matrix##
               1            2            3             4             5             6             7            8            9           10           11
1   3.735263e-16 1.393361e-16 6.777831e-17 -4.411305e-17 -1.534019e-16 -2.479452e-16  4.517006e-16 2.847301e-16 2.847301e-16 1.212307e-16 1.212307e-16
2   2.149549e-16 1.737385e-16 1.651692e-16  1.241781e-16  1.248205e-16  1.076816e-16  2.362650e-16 2.023172e-16 2.023172e-16 1.852825e-16 1.835478e-16
3   1.851877e-16 1.817280e-16 1.783694e-16  1.664480e-16  1.935578e-16  2.033205e-16  1.764348e-16 1.808873e-16 1.791526e-16 1.911944e-16 1.903270e-16
4   9.529629e-17 2.098500e-16 2.214665e-16  2.516383e-16  3.219257e-16  3.668426e-16  5.915985e-17 1.397626e-16 1.414973e-16 2.212327e-16 2.203654e-16
5   2.291139e-17 2.086793e-16 2.674814e-16  3.599316e-16  4.766682e-16  5.673842e-16 -5.412541e-17 9.972881e-17 9.625937e-17 2.483789e-16 2.483789e-16
6  -6.840130e-17 2.259194e-16 3.092574e-16  4.429735e-16  6.057463e-16  7.628802e-16 -1.541723e-16 5.030529e-17 5.030529e-17 2.808033e-16 2.773338e-16
7   5.349513e-16 1.814823e-16 9.918650e-17 -6.605555e-17 -2.534153e-16 -4.110681e-16  6.266153e-16 3.671515e-16 3.602126e-16 1.107233e-16 1.089886e-16
8   4.144261e-16 2.093539e-16 1.617721e-16  5.966987e-17 -4.243212e-17 -1.393301e-16  4.668575e-16 3.058975e-16 3.128364e-16 1.536128e-16 1.640212e-16
9   4.144261e-16 2.093539e-16 1.617721e-16  5.966987e-17 -4.243212e-17 -1.393301e-16  4.668575e-16 3.058975e-16 3.128364e-16 1.536128e-16 1.640212e-16
10  2.869612e-16 2.372251e-16 2.226231e-16  1.882141e-16  1.685511e-16  1.237329e-16  2.737060e-16 2.394410e-16 2.394410e-16 2.034404e-16 2.034404e-16
11  2.869612e-16 2.372251e-16 2.226231e-16  1.882141e-16  1.685511e-16  1.237329e-16  2.737060e-16 2.394410e-16 2.394410e-16 2.034404e-16 2.034404e-16
12  2.398458e-16 2.508530e-16 2.388591e-16  2.322187e-16  2.273131e-16  2.172047e-16  2.405472e-16 2.219125e-16 2.219125e-16 2.223631e-16 2.223631e-16
             12
1  4.975133e-17
2  1.804993e-16
3  2.018717e-16
4  2.515569e-16
5  3.006357e-16
6  3.611816e-16
7  2.402524e-17
8  1.089623e-16
9  1.089623e-16
10 1.800182e-16
11 1.800182e-16
12 2.083294e-16
##we can see H(I-H) is almost a 0 matrix##

Test for H0:β = 0 and also 95% CI for β

> sse<-t(e)%*%e
> mse<-sse/(n-k)
> se<-sqrt(mse)
> v<-solve(xtx)
> vii<-diag(v)
> t<-(b-0)/(se*sqrt(vii))
> p<-2*pt(-abs(t),n-k)
> ci<-cbind(b-se*sqrt(vii)*qt(.95,n-k),b+se*sqrt(vii)*qt(.95,n-k))
> t
        [,1]
   -5.134312
x1  6.524045
x2  7.589401
> p
           [,1]
   6.159917e-04
x1 1.083812e-04
x2 3.363518e-05
##The p-value is very small, so we reject H0:b=0 ##
> ci
##here is the CI for b ##
         [,1]       [,2]
   -64.964074 -30.780287
x1  16.470146  29.342486
x2   3.265501   5.345326

Cov[\hat\beta,e]

> covbe<-mean(b%*%t(e))-mean(b)%*%t(mean(e))
> covbe
             [,1]
[1,] 1.071988e-15
##The covariance of b and e is almost 0 ##

With Interaction Terms

Y vector and X matrix

>  dd<-read.csv("Q:\\6630a1.csv",header=T)
>  dd2<-dd
> y<-as.matrix(dd2[,4])
>  dd2$x1<-ifelse(dd2$Sex=="M",2,1)
> dd2$x2<-dd2$Educ
> dd2$x3<-dd2$x1*dd2$x2
>  xx<-as.matrix(dd2[,c(5,6,7)])
> x<-cbind(1,xx)
> dd2
    X Sex Educ  Sal x1 x2 x3
1   1   M    5 17.2  2  5 10
2   2   M    9 28.9  2  9 18
3   3   M   10 42.5  2 10 20
4   4   M   12 56.1  2 12 24
5   5   M   14 56.6  2 14 28
6   6   M   16 70.5  2 16 32
7   7   F    8 21.6  1  8  8
8   8   F   11 18.5  1 11 11
9   9   F   11 17.9  1 11 11
10 10   F   14 31.8  1 14 14
11 11   F   14 35.2  1 14 14
12 12   F   15 39.5  1 15 15
> y
## y vector ##
      [,1]
 [1,] 17.2
 [2,] 28.9
 [3,] 42.5
 [4,] 56.1
 [5,] 56.6
 [6,] 70.5
 [7,] 21.6
 [8,] 18.5
 [9,] 17.9
[10,] 31.8
[11,] 35.2
[12,] 39.5
> x
## x matrix ##
     x1 x2 x3
1  1  2  5 10
2  1  2  9 18
3  1  2 10 20
4  1  2 12 24
5  1  2 14 28
6  1  2 16 32
7  1  1  8  8
8  1  1 11 11
9  1  1 11 11
10 1  1 14 14
11 1  1 14 14
12 1  1 15 15

X'X and X'Y

>  xtx<-t(x)%*%x
> xtx
## x'x ##
        x1   x2   x3
    12  18  139  205
x1  18  30  205  337
x2 139 205 1725 2527
x3 205 337 2527 4131
> xty<- t(x)%*%y
> xty
## x'y ##
     [,1]
    436.3
x1  708.1
x2 5468.4
x3 8833.1

\hat\beta

> xtxinv<-solve(xtx)
> b<-solve(xtx,xty)
> b
## b hat ##
         [,1]
   -7.6563397
x1 -0.6527512
x2  0.9398325
x3  1.9965311

Normal Equation

> nn<-xtx%*%b
> nn
     [,1]
    436.3
x1  708.1
x2 5468.4
x3 8833.1
## So we can see the normal equation x'xb=x'y ##

residual e and \hat y

> yhat<-x%*%b
> yhat
## fitted value y hat ##
       [,1]
1  15.70263
2  35.43421
3  40.36711
4  50.23289
5  60.09868
6  69.96447
7  15.18182
8  23.99091
9  23.99091
10 32.80000
11 32.80000
12 35.73636
> e<-y-yhat
> e
## residual e ##
         [,1]
1   1.4973684
2  -6.5342105
3   2.1328947
4   5.8671053
5  -3.4986842
6   0.5355263
7   6.4181818
8  -5.4909091
9  -6.0909091
10 -1.0000000
11  2.4000000
12  3.7636364

H and (I-H)

> n<-length(y)
> k<-ncol(x)
> i<-diag(n)
> h<-x%*%xtxinv%*%t(x)
> he<-i-h
> h
## H matrix ##
               1             2             3             4             5             6             7             8             9            10            11
1   6.403509e-01  3.245614e-01  2.456140e-01  8.771930e-02 -7.017544e-02 -2.280702e-01 -2.775558e-15 -2.206568e-15 -2.206568e-15 -1.637579e-15 -1.637579e-15
2   3.245614e-01  2.192982e-01  1.929825e-01  1.403509e-01  8.771930e-02  3.508772e-02 -1.110223e-16  5.134781e-16  5.134781e-16  1.137979e-15  1.137979e-15
3   2.456140e-01  1.929825e-01  1.798246e-01  1.535088e-01  1.271930e-01  1.008772e-01 -1.998401e-15 -2.081668e-15 -2.081668e-15 -2.164935e-15 -2.164935e-15
4   8.771930e-02  1.403509e-01  1.535088e-01  1.798246e-01  2.061404e-01  2.324561e-01  4.440892e-16  6.106227e-16  6.106227e-16  7.771561e-16  7.771561e-16
5  -7.017544e-02  8.771930e-02  1.271930e-01  2.061404e-01  2.850877e-01  3.640351e-01 -1.554312e-15 -2.137179e-15 -2.137179e-15 -2.720046e-15 -2.720046e-15
6  -2.280702e-01  3.508772e-02  1.008772e-01  2.324561e-01  3.640351e-01  4.956140e-01  0.000000e+00 -6.661338e-16 -6.661338e-16 -1.332268e-15 -1.332268e-15
7  -5.329071e-15 -3.996803e-15 -3.552714e-15 -3.108624e-15 -2.664535e-15 -1.776357e-15  6.650718e-01  3.062201e-01  3.062201e-01 -5.263158e-02 -5.263158e-02
8  -3.483325e-15 -2.872702e-15 -2.747802e-15 -2.386980e-15 -2.026157e-15 -1.776357e-15  3.062201e-01  2.057416e-01  2.057416e-01  1.052632e-01  1.052632e-01
9  -3.483325e-15 -2.872702e-15 -2.747802e-15 -2.386980e-15 -2.026157e-15 -1.776357e-15  3.062201e-01  2.057416e-01  2.057416e-01  1.052632e-01  1.052632e-01
10 -1.526557e-15 -1.637579e-15 -1.720846e-15 -1.665335e-15 -1.609823e-15 -1.776357e-15 -5.263158e-02  1.052632e-01  1.052632e-01  2.631579e-01  2.631579e-01
11 -1.526557e-15 -1.637579e-15 -1.720846e-15 -1.665335e-15 -1.609823e-15 -1.776357e-15 -5.263158e-02  1.052632e-01  1.052632e-01  2.631579e-01  2.631579e-01
12 -4.996004e-16 -1.054712e-15 -1.221245e-15 -1.554312e-15 -1.887379e-15 -1.776357e-15 -1.722488e-01  7.177033e-02  7.177033e-02  3.157895e-01  3.157895e-01
              12
1  -1.373901e-15
2   1.346145e-15
3  -2.192690e-15
4   8.326673e-16
5  -2.914335e-15
6  -1.554312e-15
7  -1.722488e-01
8   7.177033e-02
9   7.177033e-02
10  3.157895e-01
11  3.157895e-01
12  3.971292e-01
> he
## I-H matrix ##
               1             2             3             4             5             6             7             8             9            10            11
1   3.596491e-01 -3.245614e-01 -2.456140e-01 -8.771930e-02  7.017544e-02  2.280702e-01  2.775558e-15  2.206568e-15  2.206568e-15  1.637579e-15  1.637579e-15
2  -3.245614e-01  7.807018e-01 -1.929825e-01 -1.403509e-01 -8.771930e-02 -3.508772e-02  1.110223e-16 -5.134781e-16 -5.134781e-16 -1.137979e-15 -1.137979e-15
3  -2.456140e-01 -1.929825e-01  8.201754e-01 -1.535088e-01 -1.271930e-01 -1.008772e-01  1.998401e-15  2.081668e-15  2.081668e-15  2.164935e-15  2.164935e-15
4  -8.771930e-02 -1.403509e-01 -1.535088e-01  8.201754e-01 -2.061404e-01 -2.324561e-01 -4.440892e-16 -6.106227e-16 -6.106227e-16 -7.771561e-16 -7.771561e-16
5   7.017544e-02 -8.771930e-02 -1.271930e-01 -2.061404e-01  7.149123e-01 -3.640351e-01  1.554312e-15  2.137179e-15  2.137179e-15  2.720046e-15  2.720046e-15
6   2.280702e-01 -3.508772e-02 -1.008772e-01 -2.324561e-01 -3.640351e-01  5.043860e-01  0.000000e+00  6.661338e-16  6.661338e-16  1.332268e-15  1.332268e-15
7   5.329071e-15  3.996803e-15  3.552714e-15  3.108624e-15  2.664535e-15  1.776357e-15  3.349282e-01 -3.062201e-01 -3.062201e-01  5.263158e-02  5.263158e-02
8   3.483325e-15  2.872702e-15  2.747802e-15  2.386980e-15  2.026157e-15  1.776357e-15 -3.062201e-01  7.942584e-01 -2.057416e-01 -1.052632e-01 -1.052632e-01
9   3.483325e-15  2.872702e-15  2.747802e-15  2.386980e-15  2.026157e-15  1.776357e-15 -3.062201e-01 -2.057416e-01  7.942584e-01 -1.052632e-01 -1.052632e-01
10  1.526557e-15  1.637579e-15  1.720846e-15  1.665335e-15  1.609823e-15  1.776357e-15  5.263158e-02 -1.052632e-01 -1.052632e-01  7.368421e-01 -2.631579e-01
11  1.526557e-15  1.637579e-15  1.720846e-15  1.665335e-15  1.609823e-15  1.776357e-15  5.263158e-02 -1.052632e-01 -1.052632e-01 -2.631579e-01  7.368421e-01
12  4.996004e-16  1.054712e-15  1.221245e-15  1.554312e-15  1.887379e-15  1.776357e-15  1.722488e-01 -7.177033e-02 -7.177033e-02 -3.157895e-01 -3.157895e-01
              12
1   1.373901e-15
2  -1.346145e-15
3   2.192690e-15
4  -8.326673e-16
5   2.914335e-15
6   1.554312e-15
7   1.722488e-01
8  -7.177033e-02
9  -7.177033e-02
10 -3.157895e-01
11 -3.157895e-01
12  6.028708e-01
> h%*%h
## Verify HH=H, H is idenpotent ##
               1             2             3             4             5             6             7             8             9            10            11
1   6.403509e-01  3.245614e-01  2.456140e-01  8.771930e-02 -7.017544e-02 -2.280702e-01 -4.944477e-15 -3.603400e-15 -3.603400e-15 -2.262323e-15 -2.262323e-15
2   3.245614e-01  2.192982e-01  1.929825e-01  1.403509e-01  8.771930e-02  3.508772e-02 -1.495879e-15 -6.169529e-16 -6.169529e-16  2.619737e-16  2.619737e-16
3   2.456140e-01  1.929825e-01  1.798246e-01  1.535088e-01  1.271930e-01  1.008772e-01 -3.190430e-15 -3.144171e-15 -3.144171e-15 -3.097912e-15 -3.097912e-15
4   8.771930e-02  1.403509e-01  1.535088e-01  1.798246e-01  2.061404e-01  2.324561e-01 -3.622833e-16 -3.160240e-16 -3.160240e-16 -2.697647e-16 -2.697647e-16
5  -7.017544e-02  8.771930e-02  1.271930e-01  2.061404e-01  2.850877e-01  3.640351e-01 -1.975028e-15 -2.927970e-15 -2.927970e-15 -3.880911e-15 -3.880911e-15
6  -2.280702e-01  3.508772e-02  1.008772e-01  2.324561e-01  3.640351e-01  4.956140e-01 -3.505967e-17 -1.321068e-15 -1.321068e-15 -2.607076e-15 -2.607076e-15
7  -1.069364e-14 -8.087541e-15 -7.368552e-15 -6.181304e-15 -4.994056e-15 -3.632572e-15  6.650718e-01  3.062201e-01  3.062201e-01 -5.263158e-02 -5.263158e-02
8  -6.922316e-15 -5.692306e-15 -5.375905e-15 -4.786664e-15 -4.197422e-15 -3.532749e-15  3.062201e-01  2.057416e-01  2.057416e-01  1.052632e-01  1.052632e-01
9  -6.922316e-15 -5.692306e-15 -5.375905e-15 -4.786664e-15 -4.197422e-15 -3.532749e-15  3.062201e-01  2.057416e-01  2.057416e-01  1.052632e-01  1.052632e-01
10 -2.973742e-15 -3.213317e-15 -3.322878e-15 -3.378389e-15 -3.433900e-15 -3.512785e-15 -5.263158e-02  1.052632e-01  1.052632e-01  2.631579e-01  2.631579e-01
11 -2.973742e-15 -3.213317e-15 -3.322878e-15 -3.378389e-15 -3.433900e-15 -3.512785e-15 -5.263158e-02  1.052632e-01  1.052632e-01  2.631579e-01  2.631579e-01
12 -1.305574e-15 -2.252186e-15 -2.558028e-15 -2.937045e-15 -3.316061e-15 -3.751385e-15 -1.722488e-01  7.177033e-02  7.177033e-02  3.157895e-01  3.157895e-01
              12
1  -1.767902e-15
2   5.789716e-16
3  -3.064313e-15
4  -2.478524e-16
5  -4.203752e-15
6  -3.052626e-15
7  -1.722488e-01
8   7.177033e-02
9   7.177033e-02
10  3.157895e-01
11  3.157895e-01
12  3.971292e-01
> he%*%he
## Verify (I-H)(I-H)=I-H, I-H is idenpotent ##
               1             2             3             4             5             6             7             8             9            10            11
1   3.596491e-01 -3.245614e-01 -2.456140e-01 -8.771930e-02  7.017544e-02  2.280702e-01  6.066386e-16  8.097368e-16  8.097368e-16  1.012835e-15  1.012835e-15
2  -3.245614e-01  7.807018e-01 -1.929825e-01 -1.403509e-01 -8.771930e-02 -3.508772e-02 -1.273835e-15 -1.643909e-15 -1.643909e-15 -2.013984e-15 -2.013984e-15
3  -2.456140e-01 -1.929825e-01  8.201754e-01 -1.535088e-01 -1.271930e-01 -1.008772e-01  8.063725e-16  1.019165e-15  1.019165e-15  1.231958e-15  1.231958e-15
4  -8.771930e-02 -1.403509e-01 -1.535088e-01  8.201754e-01 -2.061404e-01 -2.324561e-01 -1.250462e-15 -1.537269e-15 -1.537269e-15 -1.824077e-15 -1.824077e-15
5   7.017544e-02 -8.771930e-02 -1.271930e-01 -2.061404e-01  7.149123e-01 -3.640351e-01  1.133596e-15  1.346389e-15  1.346389e-15  1.559182e-15  1.559182e-15
6   2.280702e-01 -3.508772e-02 -1.008772e-01 -2.324561e-01 -3.640351e-01  5.043860e-01 -3.505967e-17  1.119962e-17  1.119962e-17  5.745891e-17  5.745891e-17
7  -3.550235e-17 -9.393514e-17 -2.631246e-16  3.594502e-17  3.350147e-16 -7.985815e-17  3.349282e-01 -3.062201e-01 -3.062201e-01  5.263158e-02  5.263158e-02
8   4.433367e-17  5.309859e-17  1.196987e-16 -1.270471e-17 -1.451081e-16  1.996454e-17 -3.062201e-01  7.942584e-01 -2.057416e-01 -1.052632e-01 -1.052632e-01
9   4.433367e-17  5.309859e-17  1.196987e-16 -1.270471e-17 -1.451081e-16  1.996454e-17 -3.062201e-01 -2.057416e-01  7.942584e-01 -1.052632e-01 -1.052632e-01
10  7.937121e-17  6.184137e-17  1.188133e-16 -4.772011e-17 -2.142536e-16  3.992907e-17  5.263158e-02 -1.052632e-01 -1.052632e-01  7.368421e-01 -2.631579e-01
11  7.937121e-17  6.184137e-17  1.188133e-16 -4.772011e-17 -2.142536e-16  3.992907e-17  5.263158e-02 -1.052632e-01 -1.052632e-01 -2.631579e-01  7.368421e-01
12 -3.063737e-16 -1.427619e-16 -1.155376e-16  1.715799e-16  4.586974e-16 -1.986715e-16  1.722488e-01 -7.177033e-02 -7.177033e-02 -3.157895e-01 -3.157895e-01
              12
1   9.799002e-16
2  -2.113319e-15
3   1.321068e-15
4  -1.913187e-15
5   1.624919e-15
6   5.599809e-17
7   1.722488e-01
8  -7.177033e-02
9  -7.177033e-02
10 -3.157895e-01
11 -3.157895e-01
12  6.028708e-01
> h%*%he
## H(I-H) is almost 0 ##
              1             2             3             4             5             6             7            8            9           10           11
1  9.854110e-16 -1.665216e-16 -4.228456e-16 -1.077741e-15 -1.684064e-15 -2.241818e-15  2.168919e-15 1.396831e-15 1.396831e-15 6.247439e-16 6.247439e-16
2  1.232640e-15  3.535829e-16  1.024311e-16 -2.905848e-16 -6.744940e-16 -1.136466e-15  1.384857e-15 1.130431e-15 1.130431e-15 8.760049e-16 8.760049e-16
3  1.269028e-15  4.726120e-16  2.770857e-16 -8.013948e-17 -4.512416e-16 -8.743853e-16  1.192029e-15 1.062503e-15 1.062503e-15 9.329769e-16 9.329769e-16
4  1.359759e-15  7.541448e-16  5.953676e-16  3.055722e-16  1.230569e-17 -3.087130e-16  8.063725e-16 9.266467e-16 9.266467e-16 1.046921e-15 1.046921e-15
5  1.455236e-15  1.022370e-15  9.063456e-16  6.742992e-16  4.630763e-16  2.518466e-16  4.207161e-16 7.907904e-16 7.907904e-16 1.160865e-15 1.160865e-15
6  1.552910e-15  1.289591e-15  1.214222e-15  1.091236e-15  9.335523e-16  8.175155e-16  3.505967e-17 6.549342e-16 6.549342e-16 1.274809e-15 1.274809e-15
7  5.364573e-15  4.090738e-15  3.815838e-15  3.072679e-15  2.329521e-15  1.856215e-15 -4.704956e-15 4.400624e-16 4.470013e-16 5.661409e-15 5.661409e-15
8  3.438991e-15  2.819603e-15  2.628103e-15  2.399684e-15  2.171265e-15  1.756392e-15  4.381041e-16 2.091473e-15 2.077595e-15 3.654635e-15 3.654635e-15
9  3.438991e-15  2.819603e-15  2.628103e-15  2.399684e-15  2.171265e-15  1.756392e-15  4.381041e-16 2.091473e-15 2.077595e-15 3.654635e-15 3.654635e-15
10 1.447185e-15  1.575738e-15  1.602032e-15  1.713055e-15  1.824077e-15  1.736428e-15  5.633721e-15 3.682242e-15 3.682242e-15 1.772400e-15 1.772400e-15
11 1.447185e-15  1.575738e-15  1.602032e-15  1.713055e-15  1.824077e-15  1.736428e-15  5.633721e-15 3.682242e-15 3.682242e-15 1.772400e-15 1.772400e-15
12 8.059741e-16  1.197474e-15  1.336783e-15  1.382732e-15  1.428682e-15  1.975028e-15  7.362105e-15 4.227492e-15 4.227492e-15 1.092876e-15 1.092876e-15
             12
1  3.940008e-16
2  7.671738e-16
3  8.716225e-16
4  1.080520e-15
5  1.289417e-15
6  1.498314e-15
7  7.302542e-15
8  4.266271e-15
9  4.266271e-15
10 1.091290e-15
11 1.091290e-15
12 7.944491e-17

Test for H0:β = 0 and also 95% CI for β

> sse<-t(e)%*%e
> mse<-sse/(n-k)
> se<-sqrt(mse)
> v<-solve(xtx)
> vii<-diag(v)
> t<-(b-0)/(se*sqrt(vii))
> p<-2*pt(-abs(t),n-k)
>  ci<-cbind(b-se*sqrt(vii)*qt(.95,n-k),b+se*sqrt(vii)*qt(.95,n-k))
> t
         [,1]
   -0.2742425
x1 -0.0414676
x2  0.4147084
x3  1.5403746
> p
        [,1]
   0.7908466
x1 0.9679394
x2 0.6892529
x3 0.1620353
## p-value is larger than 0.05, so Fail to reject H0:b=0 ## 
> ci
##The CI for b ##
          [,1]      [,2]
   -59.5714575 44.258778
x1 -29.9243295 28.618827
x2  -3.2743661  5.154031
x3  -0.4136914  4.406754

Cov[\hat\beta,e]

> covbe<-mean(b%*%t(e))-mean(b)%*%t(mean(e))
> covbe
              [,1]
[1,] -2.187118e-16
##The covariance of b&e is almost 0 ##

Jing Dong sep 22 (edt)

Proof of E[\hat\beta]=\beta, Var[\hat\beta]=\sigma^2(X'X)^{-1}

We can make the following statements without normality: E[\hat\beta] = E[(X'X)^{-1}X'Y] = (X'X)^{-1}X'E[Y] = (X'X)^{-1}X'X\beta = \beta Var[\hat\beta] = (X'X)^{-1}X'Var[Y]((X'X)^{-1}X')' = \sigma^2(X'X)^{-1}X'X(X'X)^{-1} = \sigma^2(X'X)^{-1}

Now assume that the error terms are normally distributed. Since Y = Xβ + ε is a linear combination of a multivariate normal, it is also multivariate normal. Furthermore, \hat\beta = (X'X)^{-1}X'Y is multivariate normal by the same argument. We can therefore conclude that, in the BNM, \hat\beta \sim N(\beta,\sigma^2(X'X)^{-1})

Sbush 23:57, 21 Sep 2006 (EDT)

Other

The stuff below was posted before. I suggest we organize things as above, but we can change it if you want. Can the person who was doing the last coding scheme please fill it in above?


This is the solution to the problem:

Find Cov[\hat\beta,e]

\hat\beta=(X'X)^{-1}X'Y
Var[Y] = σ2I
e = YHY = (IX(X'X) − 1X')Y
Cov[\hat\beta,e]
= (X'X) − 1X'Var[Y](IX(X'X) − 1X')'
= (X'X) − 1X2I(IX(X'X) − 1X')
= σ2(X'X) − 1X' − σ2(X'X) − 1X'X(X'X) − 1X'
= σ2(X'X) − 1X' − σ2(X'X) − 1X'
= 0

So we conclude that Cov[\hat\beta,e]=0

Jing Dong sep 23