Introduction to multivariable regression
Well, today we will look at liner regression with more than 1 predictor.  I will use example from your book (pages 241-251).

First we need to read data in (CH06FI05.DAT )

    > data<-read.table("/afs/stat.lsa.umich.edu/users/student/kutsyy/.public_html/classes/stat413.data/CH06FI05.DAT")
    > X1<-data[,1]
    > X2<-data[,2]
    > Y<-data[,3]
    > X1
     [1] 68.5 45.2 91.3 47.8 46.9 66.1 49.5 52.0 48.9 38.4 87.9 72.8 88.4 42.9 52.5
    [16] 85.7 41.3 51.7 89.6 82.7 52.3
    > X2
     [1] 16.7 16.8 18.2 16.3 17.3 18.2 15.9 17.2 16.6 16.0 18.3 17.1 17.4 15.8 17.8
    [16] 18.4 16.5 16.3 18.1 19.1 16.0
    > Y
     [1] 174.4 164.4 244.2 154.6 181.6 207.5 152.8 163.2 145.4 137.2 241.9 191.1
    [13] 232.0 145.3 161.1 209.7 146.4 144.0 232.6 224.1 166.5
Now let's create X matrix:
    > X<-cbind(1,X1,X2)
    > X
              X1   X2
     [1,] 1 68.5 16.7
     [2,] 1 45.2 16.8
     [3,] 1 91.3 18.2
     [4,] 1 47.8 16.3
     [5,] 1 46.9 17.3
     [6,] 1 66.1 18.2
     [7,] 1 49.5 15.9
     [8,] 1 52.0 17.2
     [9,] 1 48.9 16.6
    [10,] 1 38.4 16.0
    [11,] 1 87.9 18.3
    [12,] 1 72.8 17.1
    [13,] 1 88.4 17.4
    [14,] 1 42.9 15.8
    [15,] 1 52.5 17.8
    [16,] 1 85.7 18.4
    [17,] 1 41.3 16.5
    [18,] 1 51.7 16.3
    [19,] 1 89.6 18.1
    [20,] 1 82.7 19.1
    [21,] 1 52.3 16.0
Now let's try some of matrix opperation (p 242-244):
    > t(X)%*%X
                    X1       X2
         21.0  1302.40   360.00
    X1 1302.4 87707.94 22609.19
    X2  360.0 22609.19  6190.26
    > t(X)%*%Y
            [,1]
         3820.00
    X1 249643.35
    X2  66072.75
    > solve(t(X)%*%X)
                               X1           X2
        29.72892348  0.0721834719 -1.992553186
    X1   0.07218347  0.0003701761 -0.005549917
    X2  -1.99255319 -0.0055499173  0.136310637
Now let's find estimated regression function:
    > b<-solve(t(X)%*%X)%*%t(X)%*%Y
    > b
            [,1]
       -68.85707
    X1   1.45456
    X2   9.36550
Let's find fitted values and residuals:
    > cbind(Y.hat_X%*%b,e_Y-X%*%b)
              [,1]        [,2]
     [1,] 187.1841 -12.7841146
     [2,] 154.2294  10.1705737
     [3,] 234.3963   9.8036764
     [4,] 153.3285   1.2714690
     [5,] 161.3849  20.2150722
     [6,] 197.7414   9.7585779
     [7,] 152.0551   0.7449178
     [8,] 167.8666  -4.6666316
     [9,] 157.7382 -12.3381967
    [10,] 136.8460   0.3539791
    [11,] 230.3874  11.5126289
    [12,] 197.1849  -6.0849209
    [13,] 222.6857   9.3142995
    [14,] 141.5184   3.7815611
    [15,] 174.2132 -13.1132116
    [16,] 228.1239 -18.4238900
    [17,] 145.7470   0.6530062
    [18,] 159.0013 -15.0013134
    [19,] 230.9870   1.6129777
    [20,] 230.3161  -6.2160615
    [21,] 157.0644   9.4356009
Now let's try to graphicaly analyze our model:
    > motif()
    > plot(Y.hat,e)
    > plot(X1,e)
    > plot(X2,e)
    > plot(X1*X2,e)
    > plot(Y.hat,abs(e))
    > qqnorm(e)
Now let's do ANOVA:
    > t(Y)%*%Y
             [,1]
    [1,] 721072.4
    > J<-one%*%t(one)
    > one<-cbind(rep(1,21))
    > J<-one%*%t(one)
    > dim(J)
    [1] 21 21
    > SSTO<-t(Y)%*%Y-t(Y)%*%J%*%Y/21
    > SSTO
             [,1]
    [1,] 26196.21
    > SSE<-t(Y)%*%Y-t(b)%*%t(X)%*%Y
    > SSE
             [,1]
    [1,] 2180.927
    > SSR<-SSTO-SSE
    > SSR
             [,1]
    [1,] 24015.28
Test of regression and R^2:
    > F.star<-(SSR/2)/(SSE/(21-2-1))
    > F.star
            [,1]
    [1,] 99.1035
    > 1-pf(F.star,2,21-2-1)
    [1] 1.921431e-10
    > R2<-SSR/SSTO
    > R2
              [,1]
    [1,] 0.9167465
Estimating of regression parameters:
    > SB2<-solve(t(X)%*%X)*c(SSE)/18
    > SB2
                            X1           X2
        3602.03467  8.74593958 -241.4229923
    X1     8.74594  0.04485151   -0.6724426
    X2  -241.42299 -0.67244260   16.5157558
Simultaneoce Confidence intervals for beta:
    > b[1]-qt(.975,18)*sqrt(SB2[1,1])
    [1] -194.948
    > b[1]+qt(.975,18)*sqrt(SB2[1,1])
    [1] 57.23387
    > b[2]-qt(.975,18)*sqrt(SB2[2,2])
    [1] 1.009623
    > b[2]+qt(.975,18)*sqrt(SB2[2,2])
    [1] 1.899497
    > b[3]-qt(.975,18)*sqrt(SB2[3,3])
    [1] 0.8274411
    > b[3]+qt(.975,18)*sqrt(SB2[3,3])
    [1] 17.90356
Estimation of Mean Response:
    > Xh<-cbind(1,65.4,17.6)
    > Xh%*%b
             [,1]
    [1,] 191.1039
    > Xh%*%SB2%*%t(Xh)
            [,1]
    [1,] 7.65517
    > Xh%*%b-qt(.975,18)*sqrt(Xh%*%SB2%*%t(Xh))
             [,1]
    [1,] 185.2911
    > Xh%*%b+qt(.975,18)*sqrt(Xh%*%SB2%*%t(Xh))
             [,1]
    [1,] 196.9168
Prediction limits for new observasions:
    > Xh<-cbind(1,c(65.4,53.1),c(17.6,17.7))
    > Xh
         [,1] [,2] [,3]
    [1,]    1 65.4 17.6
    [2,]    1 53.1 17.7
    > 2*qf(1-.1,2,18)
    [1] 5.247894
    > qt(1-.1/4,18)
    [1] 2.100922
    > Xh%*%b
             [,1]
    [1,] 191.1039
    [2,] 174.1494
    > Xh%*%SB2%*%t(Xh)
             [,1]     [,2]
    [1,]  7.65517 10.09869
    [2,] 10.09869 21.14717
    > diag(Xh%*%SB2%*%t(Xh))
    [1]  7.65517 21.14717
    > diag(Xh%*%SB2%*%t(Xh))+SSE/18
    [1] 128.8178 142.3098
    > sqrt(diag(Xh%*%SB2%*%t(Xh))+SSE/18)
    [1] 11.34979 11.92937
    > Xh%*%b-qt(1-.1/4,18)*sqrt(diag(Xh%*%SB2%*%t(Xh))+SSE/18)
             [,1]
    [1,] 167.2589
    [2,] 149.0867
    > Xh%*%b+qt(1-.1/4,18)*sqrt(diag(Xh%*%SB2%*%t(Xh))+SSE/18)
             [,1]
    [1,] 214.9490
    [2,] 199.2121
Another way to get most of the results above:
    > fit<-lm(Y~X1+X2)
    > model.matrix(fit)
       (Intercept)   X1   X2
     1           1 68.5 16.7
     2           1 45.2 16.8
     3           1 91.3 18.2
     4           1 47.8 16.3
     5           1 46.9 17.3
     6           1 66.1 18.2
     7           1 49.5 15.9
     8           1 52.0 17.2
     9           1 48.9 16.6
    10           1 38.4 16.0
    11           1 87.9 18.3
    12           1 72.8 17.1
    13           1 88.4 17.4
    14           1 42.9 15.8
    15           1 52.5 17.8
    16           1 85.7 18.4
    17           1 41.3 16.5
    18           1 51.7 16.3
    19           1 89.6 18.1
    20           1 82.7 19.1
    21           1 52.3 16.0
    > par(mfrow=c(2,3))
    > plot(fit)
    > summary(fit)

    Call: lm(formula = Y ~ X1 + X2)
    Residuals:
        Min     1Q Median    3Q   Max
     -18.42 -6.216 0.7449 9.436 20.22

    Coefficients:
                   Value Std. Error  t value Pr(>|t|)
    (Intercept) -68.8571  60.0170    -1.1473   0.2663
             X1   1.4546   0.2118     6.8682   0.0000
             X2   9.3655   4.0640     2.3045   0.0333

    Residual standard error: 11.01 on 18 degrees of freedom
    Multiple R-Squared: 0.9167
    F-statistic: 99.1 on 2 and 18 degrees of freedom, the p-value is 1.921e-10

    Correlation of Coefficients:
       (Intercept)      X1
    X1  0.6881
    X2 -0.9898     -0.7813
    > anova(fit)
    Analysis of Variance Table

    Response: Y

    Terms added sequentially (first to last)
              Df Sum of Sq  Mean Sq  F Value      Pr(F)
    X1         1  23371.81 23371.81 192.8962 0.00000000
    X2         1    643.48   643.48   5.3108 0.03332136
    Residuals 18   2180.93   121.16
    > fit$fitted
            1        2        3        4        5        6        7        8
     187.1841 154.2294 234.3963 153.3285 161.3849 197.7414 152.0551 167.8666
            9      10       11       12       13       14       15       16
     157.7382 136.846 230.3874 197.1849 222.6857 141.5184 174.2132 228.1239
          17       18      19       20       21
     145.747 159.0013 230.987 230.3161 157.0644
    > fit$res
             1        2        3        4        5        6         7         8
     -12.78411 10.17057 9.803676 1.271469 20.21507 9.758578 0.7449178 -4.666632
            9        10       11        12       13       14        15        16
     -12.3382 0.3539791 11.51263 -6.084921 9.314299 3.781561 -13.11321 -18.42389
            17        18       19        20       21
     0.6530062 -15.00131 1.612978 -6.216062 9.435601
    > fit$coef
     (Intercept)      X1     X2
       -68.85707 1.45456 9.3655

More on working with Splus objects:
    > names(fit)
     [1] "coefficients"  "residuals"     "fitted.values" "effects"
     [5] "R"             "rank"          "assign"        "df.residual"
     [9] "contrasts"     "terms"         "call"
    >  names(summary(fit))
     [1] "call"         "terms"        "residuals"    "coefficients" "sigma"
     [6] "df"           "r.squared"    "fstatistic"   "cov.unscaled" "correlation"
    > names(anova(fit)) names(anova(fit))
    [1] "Df"        "Sum of Sq" "Mean Sq"   "F Value"   "Pr(F)"
    > anova(fit)$Pr
    [1] 4.639555e-11 3.332136e-02           NA
    > summary(fit)$sigma
    [1] 11.00739