R version 3.6.3 (2020-02-29) -- "Holding the Windsock"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> pdf("g04.pdf");options(width=64)
> #Examine the parent-child height data.  Data are from
> #https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/T0HSJ1
> #Documentation claims that these are formatted for Stata; fortunately R reads 
> #them just fine.  Tell R that the first row contains column names.
> galton<-read.table("galton-stata11.tab",header=TRUE)
> #First examine daughters.
> daughters<-galton[galton$female==1,]
> #Families have multiple children.  Most of what we will do this semester will
> #require observations to be independent, and so keep only the first daughter.
> first<-c(TRUE,daughters$family[-1]!=daughters$family[-length(daughters$family)])
> fd<-daughters[first,c("mother","father","height")]
> fd$mothercm<-fd$mother*2.54
> attach(fd)
> Z<-cbind(1,mother-mean(mother))
> nextz<-fd$mothercm-Z%*%solve(t(Z)%*%Z)%*%t(Z)%*%fd$mothercm
> nextz
                [,1]
  [1,] -1.136868e-13
  [2,] -2.273737e-13
  [3,] -1.136868e-13
  [4,] -1.136868e-13
  [5,]  2.842171e-14
  [6,] -3.410605e-13
  [7,] -3.410605e-13
  [8,] -2.273737e-13
  [9,]  2.273737e-13
 [10,] -2.842171e-14
 [11,]  1.136868e-13
 [12,] -8.526513e-14
 [13,] -2.273737e-13
 [14,] -1.705303e-13
 [15,]  2.842171e-14
 [16,] -1.136868e-13
 [17,] -5.684342e-14
 [18,]  1.136868e-13
 [19,] -3.410605e-13
 [20,] -1.136868e-13
 [21,] -1.705303e-13
 [22,] -2.842171e-14
 [23,] -1.136868e-13
 [24,] -5.684342e-14
 [25,] -5.684342e-14
 [26,] -5.684342e-14
 [27,]  1.136868e-13
 [28,]  1.136868e-13
 [29,]  1.136868e-13
 [30,]  1.136868e-13
 [31,]  1.136868e-13
 [32,] -8.526513e-14
 [33,]  1.136868e-13
 [34,] -1.136868e-13
 [35,]  2.273737e-13
 [36,]  2.273737e-13
 [37,]  2.273737e-13
 [38,]  2.273737e-13
 [39,] -2.842171e-14
 [40,] -2.842171e-14
 [41,] -1.705303e-13
 [42,] -1.136868e-13
 [43,]  2.842171e-14
 [44,]  2.842171e-14
 [45,] -5.684342e-14
 [46,]  1.136868e-13
 [47,] -5.684342e-14
 [48,] -5.684342e-14
 [49,]  1.136868e-13
 [50,]  1.136868e-13
 [51,] -2.273737e-13
 [52,]  1.136868e-13
 [53,] -8.526513e-14
 [54,]  1.136868e-13
 [55,]  1.136868e-13
 [56,] -1.136868e-13
 [57,] -1.136868e-13
 [58,] -2.273737e-13
 [59,] -1.705303e-13
 [60,] -1.705303e-13
 [61,] -1.705303e-13
 [62,] -1.705303e-13
 [63,] -1.705303e-13
 [64,] -1.705303e-13
 [65,]  2.842171e-14
 [66,] -1.136868e-13
 [67,] -1.136868e-13
 [68,]  0.000000e+00
 [69,] -1.136868e-13
 [70,]  2.842171e-14
 [71,] -1.136868e-13
 [72,] -1.136868e-13
 [73,] -8.526513e-14
 [74,] -5.684342e-14
 [75,]  1.136868e-13
 [76,] -5.684342e-14
 [77,] -5.684342e-14
 [78,]  1.136868e-13
 [79,]  5.684342e-14
 [80,]  1.136868e-13
 [81,]  2.842171e-14
 [82,]  2.842171e-14
 [83,]  2.842171e-14
 [84,]  2.842171e-14
 [85,] -5.684342e-14
 [86,] -1.136868e-13
 [87,]  2.273737e-13
 [88,]  2.273737e-13
 [89,] -2.273737e-13
 [90,] -2.273737e-13
 [91,] -2.273737e-13
 [92,]  2.273737e-13
 [93,]  2.273737e-13
 [94,] -1.705303e-13
 [95,]  2.842171e-14
 [96,] -1.136868e-13
 [97,] -5.684342e-14
 [98,] -5.684342e-14
 [99,]  1.136868e-13
[100,]  1.136868e-13
[101,]  1.136868e-13
[102,]  1.136868e-13
[103,] -2.273737e-13
[104,]  1.136868e-13
[105,] -8.526513e-14
[106,] -8.526513e-14
[107,]  2.842171e-14
[108,]  2.842171e-14
[109,] -5.684342e-14
[110,] -1.136868e-13
[111,] -2.273737e-13
[112,] -2.842171e-14
[113,] -2.842171e-14
[114,] -1.705303e-13
[115,] -1.705303e-13
[116,] -1.136868e-13
[117,] -1.136868e-13
[118,] -1.136868e-13
[119,]  2.842171e-14
[120,] -1.136868e-13
[121,] -5.684342e-14
[122,]  1.136868e-13
[123,] -5.684342e-14
[124,] -5.684342e-14
[125,] -5.684342e-14
[126,]  1.136868e-13
[127,]  1.136868e-13
[128,] -8.526513e-14
[129,]  2.842171e-14
[130,]  2.842171e-14
[131,]  2.842171e-14
[132,]  1.421085e-13
[133,] -1.705303e-13
[134,] -2.842171e-14
[135,] -2.842171e-14
[136,] -1.705303e-13
[137,] -1.705303e-13
[138,]  1.136868e-13
[139,] -5.684342e-14
[140,]  1.136868e-13
[141,] -1.136868e-13
[142,] -1.136868e-13
[143,]  2.273737e-13
[144,]  2.273737e-13
[145,] -1.705303e-13
[146,] -2.842171e-14
[147,]  1.136868e-13
[148,] -5.684342e-14
[149,] -2.273737e-13
[150,]  2.842171e-14
[151,]  2.842171e-14
[152,] -1.136868e-13
[153,] -1.136868e-13
[154,]  2.273737e-13
[155,] -1.705303e-13
[156,] -2.842171e-14
[157,] -1.705303e-13
[158,] -1.136868e-13
[159,] -5.684342e-14
[160,] -5.684342e-14
[161,] -5.684342e-14
[162,]  2.842171e-14
[163,] -1.136868e-13
[164,] -1.136868e-13
[165,]  2.842171e-14
[166,]  2.273737e-13
[167,] -5.684342e-14
[168,] -1.705303e-13
> X<-cbind(1,fd$mother,fd$mothercm)
> M<-t(X)%*%X
> M
         [,1]      [,2]       [,3]
[1,]   168.00   10770.5   27357.07
[2,] 10770.50  691393.1 1756138.60
[3,] 27357.07 1756138.6 4460592.05
> solve(M)
             [,1]          [,2]          [,3]
[1,]  4.597672414 -4.857063e-02 -9.074437e-03
[2,] -0.050320181 -5.109703e+08  2.011694e+08
[3,] -0.008386697  2.011694e+08 -7.920056e+07
> summary(lm(height~mother+mothercm,data=fd))

Call:
lm(formula = height ~ mother + mothercm, data = fd)

Residuals:
   Min     1Q Median     3Q    Max 
-7.708 -1.221  0.096  1.292  5.792 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 42.24854    4.67281   9.041 3.92e-16 ***
mother       0.35651    0.07284   4.894 2.32e-06 ***
mothercm          NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.179 on 166 degrees of freedom
Multiple R-squared:  0.1261,	Adjusted R-squared:  0.1208 
F-statistic: 23.95 on 1 and 166 DF,  p-value: 2.32e-06

> fd$mothercmround<-round(fd$mother*2.54)
> summary(lm(height~mother,data=fd))

Call:
lm(formula = height ~ mother, data = fd)

Residuals:
   Min     1Q Median     3Q    Max 
-7.708 -1.221  0.096  1.292  5.792 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 42.24854    4.67281   9.041 3.92e-16 ***
mother       0.35651    0.07284   4.894 2.32e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.179 on 166 degrees of freedom
Multiple R-squared:  0.1261,	Adjusted R-squared:  0.1208 
F-statistic: 23.95 on 1 and 166 DF,  p-value: 2.32e-06

> summary(colinfit<-lm(height~mother+mothercmround,data=fd))

Call:
lm(formula = height ~ mother + mothercmround, data = fd)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.7377 -1.2562  0.1006  1.3926  5.7623 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    43.3121     4.7518   9.115 2.59e-16 ***
mother         -1.4345     1.5073  -0.952    0.343    
mothercmround   0.6988     0.5874   1.190    0.236    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.177 on 165 degrees of freedom
Multiple R-squared:  0.1335,	Adjusted R-squared:  0.123 
F-statistic: 12.71 on 2 and 165 DF,  p-value: 7.316e-06

> # First install via install.packages("regclass")
> library(regclass)# For VIF
> VIF(colinfit)
       mother mothercmround 
     429.3064      429.3064 
> # Exhibit residual standard deviation fit.
> summary(myfit<-lm(height~mother+father,data=fd))

Call:
lm(formula = height ~ mother + father, data = fd)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6346 -1.0986 -0.0422  1.3220  4.7772 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 17.10803    5.43661   3.147  0.00196 ** 
mother       0.32877    0.06415   5.125 8.24e-07 ***
father       0.38753    0.05490   7.059 4.43e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.916 on 165 degrees of freedom
Multiple R-squared:  0.3288,	Adjusted R-squared:  0.3207 
F-statistic: 40.41 on 2 and 165 DF,  p-value: 5.19e-15

> VIF(myfit)
  mother   father 
1.003766 1.003766 
> X<-cbind(1,fd$mother,fd$father)
> H<-diag(rep(1,dim(X)[1]))-X%*%solve(t(X)%*%X)%*%t(X)
> fd$height%*%H%*%fd$height/(dim(X)[1]-dim(X)[2])
         [,1]
[1,] 3.669732
> X<-cbind(1,fd$mother,fd$father)
> est<-solve(t(X)%*%X)%*%t(X)%*%fd$height
> varmat<-solve(t(X)%*%X)
> est[2:3]%*%solve(varmat[2:3,2:3])%*%est[2:3]
         [,1]
[1,] 296.6213
> # This R code doesn't give you quite what you might expect.
> anova(myfit)
Analysis of Variance Table

Response: height
           Df Sum Sq Mean Sq F value    Pr(>F)    
mother      1 113.77  113.77  31.001 1.026e-07 ***
father      1 182.86  182.86  49.828 4.431e-11 ***
Residuals 165 605.51    3.67                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> # Library rms must have been previously installed via 
> # install.packages("rms")
> library(rms)# For ols
> anova(ols(height~mother+father,data=fd))
                Analysis of Variance          Response: height 

 Factor     d.f. Partial SS MS         F     P     
 mother       1   96.38921   96.389208 26.27 <.0001
 father       1  182.85570  182.855696 49.83 <.0001
 REGRESSION   2  296.62131  148.310653 40.41 <.0001
 ERROR      165  605.50578    3.669732             
> # This gives something closer to what is traditionally produced, except that
> # there's no total line, and each of the separate explanatory variables is
> # given.
> library(NonparametricHeuristic)# For traditionalanova
> traditionalanova(myfit)
         Sum of Squares Degrees of Freedom Mean Squares
Model          296.6213                  2   148.310653
Residual       605.5058                165     3.669732
Total          902.1271                167           NA
                F            P
Model    40.41457 5.190056e-15
Residual       NA           NA
Total          NA           NA
> 
<\/pre>