src/test/R/multipleOLSRegressionTestCases - commons-math - Git at Google

 # Licensed to the Apache Software Foundation (ASF) under one or more
 # contributor license agreements.  See the NOTICE file distributed with
 # this work for additional information regarding copyright ownership.
 # The ASF licenses this file to You under the Apache License, Version 2.0
 # (the "License"); you may not use this file except in compliance with
 # the License.  You may obtain a copy of the License at
 #
 #      http://www.apache.org/licenses/LICENSE-2.0
 #
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
 #
 #------------------------------------------------------------------------------
 # R source file to validate OLS multiple regression tests in
 # org.apache.commons.math.stat.regression.OLSMultipleLinearRegressionTest
 #
 # To run the test, install R, put this file and testFunctions
 # into the same directory, launch R from this directory and then enter
 # source("<name-of-this-file>")
 #
 #------------------------------------------------------------------------------
 tol <- 1E-8                       # default error tolerance for tests
 #------------------------------------------------------------------------------
 # Function definitions

 source("testFunctions")           # utility test functions
 options(digits=16)	              # override number of digits displayed

 # function to verify OLS computations

 verifyRegression <- function(model, expectedBeta, expectedResiduals,
   expectedErrors, expectedStdError, expectedRSquare, expecteAdjRSquare, modelName) {
     betaHat <- as.vector(coefficients(model))
     residuals <- as.vector(residuals(model))
     errors <-  as.vector(as.matrix(coefficients(summary(model)))[,2])
     stdError <- summary(model)$sigma
     rSquare <- summary(model)$r.squared
     adjRSquare <- summary(model)$adj.r.squared
     output <- c("Parameter test dataset = ", modelName)
     if (assertEquals(expectedBeta,betaHat,tol,"Parameters")) {
         displayPadded(output, SUCCEEDED, WIDTH)
     } else {
         displayPadded(output, FAILED, WIDTH)
     }
     output <- c("Residuals test dataset = ", modelName)
     if (assertEquals(expectedResiduals,residuals,tol,"Residuals")) {
         displayPadded(output, SUCCEEDED, WIDTH)
     } else {
         displayPadded(output, FAILED, WIDTH)
     }
     output <- c("Errors test dataset = ", modelName)
     if (assertEquals(expectedErrors,errors,tol,"Errors")) {
         displayPadded(output, SUCCEEDED, WIDTH)
     } else {
         displayPadded(output, FAILED, WIDTH)
     }
     output <- c("Standard Error test dataset = ", modelName)
     if (assertEquals(expectedStdError,stdError,tol,"Regression Standard Error")) {
         displayPadded(output, SUCCEEDED, WIDTH)
     } else {
         displayPadded(output, FAILED, WIDTH)
     }
     output <- c("RSquared test dataset = ", modelName)
     if (assertEquals(expectedRSquare,rSquare,tol,"RSquared")) {
         displayPadded(output, SUCCEEDED, WIDTH)
     } else {
         displayPadded(output, FAILED, WIDTH)
     }
     output <- c("Adjusted RSquared test dataset = ", modelName)
     if (assertEquals(expecteAdjRSquare,adjRSquare,tol,"Adjusted RSquared")) {
         displayPadded(output, SUCCEEDED, WIDTH)
     } else {
         displayPadded(output, FAILED, WIDTH)
     }
 }

 #--------------------------------------------------------------------------
 cat("Multiple regression OLS test cases\n")

 # Perfect fit
 x1 <- c(0,2,0,0,0,0)
 x2 <- c(0,0,3,0,0,0)
 x3 <- c(0,0,0,4,0,0)
 x4 <- c(0,0,0,0,5,0)
 x5 <- c(0,0,0,0,0,6)
 y <- c(11, 12, 13, 14, 15, 16)
 model <- lm(y ~ x1 + x2 + x3 + x4 + x5)
 expectedBeta <- c(11.0,0.5,0.666666666666667,0.75,0.8,0.8333333333333333)
 expectedResiduals <- c(0,0,0,0,0,0)
 expectedErrors <- c(NaN,NaN,NaN,NaN,NaN,NaN)
 expectedStdError <- NaN
 expectedRSquare <- 1
 expectedAdjRSquare <- NaN
 verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
  expectedStdError, expectedRSquare, expectedAdjRSquare, "perfect fit")

 # Longley
 #
 # Data Source: J. Longley (1967) "An Appraisal of Least Squares Programs for the
 # Electronic Computer from the Point of View of the User",
 # Journal of the American Statistical Association,
 # vol. 62. September, pp. 819-841.
 #
 # Certified values (and data) are from NIST:
 # http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Longley.dat
 #
 design <- matrix(c(60323,83.0,234289,2356,1590,107608,1947,
                     61122,88.5,259426,2325,1456,108632,1948,
                     60171,88.2,258054,3682,1616,109773,1949,
                     61187,89.5,284599,3351,1650,110929,1950,
                     63221,96.2,328975,2099,3099,112075,1951,
                     63639,98.1,346999,1932,3594,113270,1952,
                     64989,99.0,365385,1870,3547,115094,1953,
                     63761,100.0,363112,3578,3350,116219,1954,
                     66019,101.2,397469,2904,3048,117388,1955,
                     67857,104.6,419180,2822,2857,118734,1956,
                     68169,108.4,442769,2936,2798,120445,1957,
                     66513,110.8,444546,4681,2637,121950,1958,
                     68655,112.6,482704,3813,2552,123366,1959,
                     69564,114.2,502601,3931,2514,125368,1960,
                     69331,115.7,518173,4806,2572,127852,1961,
                     70551,116.9,554894,4007,2827,130081,1962),
                     nrow = 16, ncol = 7, byrow = TRUE)
 y <- design[,1]
 x1 <- design[,2]
 x2 <- design[,3]
 x3 <- design[,4]
 x4 <- design[,5]
 x5 <- design[,6]
 x6 <- design[,7]
 model <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6)

 estimates <- matrix(c(-3482258.63459582,890420.383607373,
                        15.0618722713733,84.9149257747669,
                       -0.358191792925910E-01,0.334910077722432E-01,
                       -2.02022980381683,0.488399681651699,
                       -1.03322686717359,0.214274163161675,
                       -0.511041056535807E-01,0.226073200069370,
                        1829.15146461355,455.478499142212),
                        nrow = 7, ncol = 2, byrow = TRUE)

 expectedBeta <- estimates[,1]
 expectedErrors <- estimates[,2]
 expectedResiduals <- c(267.340029759711,-94.0139423988359,46.28716775752924,
  -410.114621930906,309.7145907602313,-249.3112153297231,-164.0489563956039,
  -13.18035686637081,14.30477260005235,455.394094551857,-17.26892711483297,
  -39.0550425226967,-155.5499735953195,-85.6713080421283,341.9315139607727,
  -206.7578251937366)
 expectedStdError <- 304.8540735619638
 expectedRSquare <- 0.995479004577296
 expectedAdjRSquare <- 0.992465007628826

 verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
 expectedStdError, expectedRSquare, expectedAdjRSquare, "Longley")

 # Model with no intercept
 model <- lm(y ~ 0 + x1 + x2 + x3 + x4 + x5 + x6)

 estimates <- matrix(c(-52.99357013868291, 129.54486693117232,
          0.07107319907358, 0.03016640003786,
         -0.42346585566399, 0.41773654056612,
         -0.57256866841929, 0.27899087467676,
         -0.41420358884978, 0.32128496193363,
          48.41786562001326, 17.68948737819961),
          nrow = 6, ncol = 2, byrow = TRUE)

 expectedBeta <- estimates[,1]
 expectedErrors <- estimates[,2]
 expectedResiduals <- c(279.90274927293092, -130.32465380836874, 90.73228661967445,
   -401.31252201634948, -440.46768772620027, -543.54512853774793, 201.32111639536299,
    215.90889365977932, 73.09368242049943, 913.21694494481869, 424.82484953610174,
    -8.56475876776709, -361.32974610842876, 27.34560497213464, 151.28955976355002,
    -492.49937355336846)
 expectedStdError <- 475.1655079819517
 expectedRSquare <- 0.9999670130706
 expectedAdjRSquare <- 0.999947220913

 verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
 expectedStdError, expectedRSquare, expectedAdjRSquare, "Longley No Intercept")

 # Swiss Fertility (R dataset named "swiss")

 design <- matrix(c(80.2,17.0,15,12,9.96,
   83.1,45.1,6,9,84.84,
   92.5,39.7,5,5,93.40,
   85.8,36.5,12,7,33.77,
   76.9,43.5,17,15,5.16,
   76.1,35.3,9,7,90.57,
   83.8,70.2,16,7,92.85,
   92.4,67.8,14,8,97.16,
   82.4,53.3,12,7,97.67,
   82.9,45.2,16,13,91.38,
   87.1,64.5,14,6,98.61,
   64.1,62.0,21,12,8.52,
   66.9,67.5,14,7,2.27,
   68.9,60.7,19,12,4.43,
   61.7,69.3,22,5,2.82,
   68.3,72.6,18,2,24.20,
   71.7,34.0,17,8,3.30,
   55.7,19.4,26,28,12.11,
   54.3,15.2,31,20,2.15,
   65.1,73.0,19,9,2.84,
   65.5,59.8,22,10,5.23,
   65.0,55.1,14,3,4.52,
   56.6,50.9,22,12,15.14,
   57.4,54.1,20,6,4.20,
   72.5,71.2,12,1,2.40,
   74.2,58.1,14,8,5.23,
   72.0,63.5,6,3,2.56,
   60.5,60.8,16,10,7.72,
   58.3,26.8,25,19,18.46,
   65.4,49.5,15,8,6.10,
   75.5,85.9,3,2,99.71,
   69.3,84.9,7,6,99.68,
   77.3,89.7,5,2,100.00,
   70.5,78.2,12,6,98.96,
   79.4,64.9,7,3,98.22,
   65.0,75.9,9,9,99.06,
   92.2,84.6,3,3,99.46,
   79.3,63.1,13,13,96.83,
   70.4,38.4,26,12,5.62,
   65.7,7.7,29,11,13.79,
   72.7,16.7,22,13,11.22,
   64.4,17.6,35,32,16.92,
   77.6,37.6,15,7,4.97,
   67.6,18.7,25,7,8.65,
   35.0,1.2,37,53,42.34,
   44.7,46.6,16,29,50.43,
   42.8,27.7,22,29,58.33),
   nrow = 47, ncol = 5, byrow = TRUE)

 y  <- design[,1]
 x1 <- design[,2]
 x2 <- design[,3]
 x3 <- design[,4]
 x4 <- design[,5]

 model <- lm(y ~ x1 + x2 + x3 + x4)

 estimates <- matrix(c(91.05542390271397,6.94881329475087,
                       -0.22064551045715,0.07360008972340,
                       -0.26058239824328,0.27410957467466,
                       -0.96161238456030,0.19454551679325,
                        0.12441843147162,0.03726654773803),
                        nrow = 5, ncol = 2, byrow = TRUE)

 expectedBeta <- estimates[,1]
 expectedErrors <- estimates[,2]

 expectedResiduals <- c(7.1044267859730512,1.6580347433531366,
   4.6944952770029644,8.4548022690166160,13.6547432343186212,
  -9.3586864458500774,7.5822446330520386,15.5568995563859289,
   0.8113090736598980,7.1186762732484308,7.4251378771228724,
   2.6761316873234109,0.8351584810309354,7.1769991119615177,
  -3.8746753206299553,-3.1337779476387251,-0.1412575244091504,
   1.1186809170469780,-6.3588097346816594,3.4039270429434074,
   2.3374058329820175,-7.9272368576900503,-7.8361010968497959,
  -11.2597369269357070,0.9445333697827101,6.6544245101380328,
  -0.9146136301118665,-4.3152449403848570,-4.3536932047009183,
  -3.8907885169304661,-6.3027643926302188,-7.8308982189289091,
  -3.1792280015332750,-6.7167298771158226,-4.8469946718041754,
  -10.6335664353633685,11.1031134362036958,6.0084032641811733,
   5.4326230830188482,-7.2375578629692230,2.1671550814448222,
   15.0147574652763112,4.8625103516321015,-7.1597256413907706,
   -0.4515205619767598,-10.2916870903837587,-15.7812984571900063)

 expectedStdError <- 7.73642194433223
 expectedRSquare <- 0.649789742860228
 expectedAdjRSquare <- 0.6164363850373927

 verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
 expectedStdError, expectedRSquare, expectedAdjRSquare, "Swiss Fertility")

 # model with no intercept
 model <- lm(y ~ 0 + x1 + x2 + x3 + x4)

 estimates <- matrix(c(0.52191832900513, 0.10470063765677,
       2.36588087917963, 0.41684100584290,
      -0.94770353802795, 0.43370143099691,
       0.30851985863609, 0.07694953606522),
      nrow = 4, ncol = 2, byrow = TRUE)

 expectedBeta <- estimates[,1]
 expectedErrors <- estimates[,2]

 expectedResiduals <- c(44.138759883538249, 27.720705122356215, 35.873200836126799,
   34.574619581211977, 26.600168342080213, 15.074636243026923, -12.704904871199814,
   1.497443824078134, 2.691972687079431, 5.582798774291231, -4.422986561283165,
   -9.198581600334345, 4.481765170730647, 2.273520207553216, -22.649827853221336,
   -17.747900013943308, 20.298314638496436, 6.861405135329779, -8.684712790954924,
   -10.298639278062371, -9.896618896845819, 4.568568616351242, -15.313570491727944,
   -13.762961360873966, 7.156100301980509, 16.722282219843990, 26.716200609071898,
   -1.991466398777079, -2.523342564719335, 9.776486693095093, -5.297535127628603,
   -16.639070567471094, -10.302057295211819, -23.549487860816846, 1.506624392156384,
   -17.939174438345930, 13.105792202765040, -1.943329906928462, -1.516005841666695,
   -0.759066561832886, 20.793137744128977, -2.485236153005426, 27.588238710486976,
   2.658333257106881, -15.998337823623046, -5.550742066720694, -14.219077806826615)

 expectedStdError <- 17.24710630547
 expectedRSquare <- 0.946350722085
 expectedAdjRSquare <- 0.9413600915813

 verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
 expectedStdError, expectedRSquare, expectedAdjRSquare, "Swiss Fertility No Intercept")

 displayDashes(WIDTH)
	# Licensed to the Apache Software Foundation (ASF) under one or more
	# contributor license agreements. See the NOTICE file distributed with
	# this work for additional information regarding copyright ownership.
	# The ASF licenses this file to You under the Apache License, Version 2.0
	# (the "License"); you may not use this file except in compliance with
	# the License. You may obtain a copy of the License at
	#
	# http://www.apache.org/licenses/LICENSE-2.0
	#
	# Unless required by applicable law or agreed to in writing, software
	# distributed under the License is distributed on an "AS IS" BASIS,
	# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	# See the License for the specific language governing permissions and
	# limitations under the License.
	#
	#------------------------------------------------------------------------------
	# R source file to validate OLS multiple regression tests in
	# org.apache.commons.math.stat.regression.OLSMultipleLinearRegressionTest
	#
	# To run the test, install R, put this file and testFunctions
	# into the same directory, launch R from this directory and then enter
	# source("<name-of-this-file>")
	#
	#------------------------------------------------------------------------------
	tol <- 1E-8 # default error tolerance for tests
	#------------------------------------------------------------------------------
	# Function definitions

	source("testFunctions") # utility test functions
	options(digits=16) # override number of digits displayed

	# function to verify OLS computations

	verifyRegression <- function(model, expectedBeta, expectedResiduals,
	expectedErrors, expectedStdError, expectedRSquare, expecteAdjRSquare, modelName) {
	betaHat <- as.vector(coefficients(model))
	residuals <- as.vector(residuals(model))
	errors <- as.vector(as.matrix(coefficients(summary(model)))[,2])
	stdError <- summary(model)$sigma
	rSquare <- summary(model)$r.squared
	adjRSquare <- summary(model)$adj.r.squared
	output <- c("Parameter test dataset = ", modelName)
	if (assertEquals(expectedBeta,betaHat,tol,"Parameters")) {
	displayPadded(output, SUCCEEDED, WIDTH)
	} else {
	displayPadded(output, FAILED, WIDTH)
	}
	output <- c("Residuals test dataset = ", modelName)
	if (assertEquals(expectedResiduals,residuals,tol,"Residuals")) {
	displayPadded(output, SUCCEEDED, WIDTH)
	} else {
	displayPadded(output, FAILED, WIDTH)
	}
	output <- c("Errors test dataset = ", modelName)
	if (assertEquals(expectedErrors,errors,tol,"Errors")) {
	displayPadded(output, SUCCEEDED, WIDTH)
	} else {
	displayPadded(output, FAILED, WIDTH)
	}
	output <- c("Standard Error test dataset = ", modelName)
	if (assertEquals(expectedStdError,stdError,tol,"Regression Standard Error")) {
	displayPadded(output, SUCCEEDED, WIDTH)
	} else {
	displayPadded(output, FAILED, WIDTH)
	}
	output <- c("RSquared test dataset = ", modelName)
	if (assertEquals(expectedRSquare,rSquare,tol,"RSquared")) {
	displayPadded(output, SUCCEEDED, WIDTH)
	} else {
	displayPadded(output, FAILED, WIDTH)
	}
	output <- c("Adjusted RSquared test dataset = ", modelName)
	if (assertEquals(expecteAdjRSquare,adjRSquare,tol,"Adjusted RSquared")) {
	displayPadded(output, SUCCEEDED, WIDTH)
	} else {
	displayPadded(output, FAILED, WIDTH)
	}
	}

	#--------------------------------------------------------------------------
	cat("Multiple regression OLS test cases\n")

	# Perfect fit
	x1 <- c(0,2,0,0,0,0)
	x2 <- c(0,0,3,0,0,0)
	x3 <- c(0,0,0,4,0,0)
	x4 <- c(0,0,0,0,5,0)
	x5 <- c(0,0,0,0,0,6)
	y <- c(11, 12, 13, 14, 15, 16)
	model <- lm(y ~ x1 + x2 + x3 + x4 + x5)
	expectedBeta <- c(11.0,0.5,0.666666666666667,0.75,0.8,0.8333333333333333)
	expectedResiduals <- c(0,0,0,0,0,0)
	expectedErrors <- c(NaN,NaN,NaN,NaN,NaN,NaN)
	expectedStdError <- NaN
	expectedRSquare <- 1
	expectedAdjRSquare <- NaN
	verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
	expectedStdError, expectedRSquare, expectedAdjRSquare, "perfect fit")

	# Longley
	#
	# Data Source: J. Longley (1967) "An Appraisal of Least Squares Programs for the
	# Electronic Computer from the Point of View of the User",
	# Journal of the American Statistical Association,
	# vol. 62. September, pp. 819-841.
	#
	# Certified values (and data) are from NIST:
	# http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Longley.dat
	#
	design <- matrix(c(60323,83.0,234289,2356,1590,107608,1947,
	61122,88.5,259426,2325,1456,108632,1948,
	60171,88.2,258054,3682,1616,109773,1949,
	61187,89.5,284599,3351,1650,110929,1950,
	63221,96.2,328975,2099,3099,112075,1951,
	63639,98.1,346999,1932,3594,113270,1952,
	64989,99.0,365385,1870,3547,115094,1953,
	63761,100.0,363112,3578,3350,116219,1954,
	66019,101.2,397469,2904,3048,117388,1955,
	67857,104.6,419180,2822,2857,118734,1956,
	68169,108.4,442769,2936,2798,120445,1957,
	66513,110.8,444546,4681,2637,121950,1958,
	68655,112.6,482704,3813,2552,123366,1959,
	69564,114.2,502601,3931,2514,125368,1960,
	69331,115.7,518173,4806,2572,127852,1961,
	70551,116.9,554894,4007,2827,130081,1962),
	nrow = 16, ncol = 7, byrow = TRUE)
	y <- design[,1]
	x1 <- design[,2]
	x2 <- design[,3]
	x3 <- design[,4]
	x4 <- design[,5]
	x5 <- design[,6]
	x6 <- design[,7]
	model <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6)

	estimates <- matrix(c(-3482258.63459582,890420.383607373,
	15.0618722713733,84.9149257747669,
	-0.358191792925910E-01,0.334910077722432E-01,
	-2.02022980381683,0.488399681651699,
	-1.03322686717359,0.214274163161675,
	-0.511041056535807E-01,0.226073200069370,
	1829.15146461355,455.478499142212),
	nrow = 7, ncol = 2, byrow = TRUE)

	expectedBeta <- estimates[,1]
	expectedErrors <- estimates[,2]
	expectedResiduals <- c(267.340029759711,-94.0139423988359,46.28716775752924,
	-410.114621930906,309.7145907602313,-249.3112153297231,-164.0489563956039,
	-13.18035686637081,14.30477260005235,455.394094551857,-17.26892711483297,
	-39.0550425226967,-155.5499735953195,-85.6713080421283,341.9315139607727,
	-206.7578251937366)
	expectedStdError <- 304.8540735619638
	expectedRSquare <- 0.995479004577296
	expectedAdjRSquare <- 0.992465007628826

	verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
	expectedStdError, expectedRSquare, expectedAdjRSquare, "Longley")

	# Model with no intercept
	model <- lm(y ~ 0 + x1 + x2 + x3 + x4 + x5 + x6)

	estimates <- matrix(c(-52.99357013868291, 129.54486693117232,
	0.07107319907358, 0.03016640003786,
	-0.42346585566399, 0.41773654056612,
	-0.57256866841929, 0.27899087467676,
	-0.41420358884978, 0.32128496193363,
	48.41786562001326, 17.68948737819961),
	nrow = 6, ncol = 2, byrow = TRUE)

	expectedBeta <- estimates[,1]
	expectedErrors <- estimates[,2]
	expectedResiduals <- c(279.90274927293092, -130.32465380836874, 90.73228661967445,
	-401.31252201634948, -440.46768772620027, -543.54512853774793, 201.32111639536299,
	215.90889365977932, 73.09368242049943, 913.21694494481869, 424.82484953610174,
	-8.56475876776709, -361.32974610842876, 27.34560497213464, 151.28955976355002,
	-492.49937355336846)
	expectedStdError <- 475.1655079819517
	expectedRSquare <- 0.9999670130706
	expectedAdjRSquare <- 0.999947220913

	verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
	expectedStdError, expectedRSquare, expectedAdjRSquare, "Longley No Intercept")

	# Swiss Fertility (R dataset named "swiss")

	design <- matrix(c(80.2,17.0,15,12,9.96,
	83.1,45.1,6,9,84.84,
	92.5,39.7,5,5,93.40,
	85.8,36.5,12,7,33.77,
	76.9,43.5,17,15,5.16,
	76.1,35.3,9,7,90.57,
	83.8,70.2,16,7,92.85,
	92.4,67.8,14,8,97.16,
	82.4,53.3,12,7,97.67,
	82.9,45.2,16,13,91.38,
	87.1,64.5,14,6,98.61,
	64.1,62.0,21,12,8.52,
	66.9,67.5,14,7,2.27,
	68.9,60.7,19,12,4.43,
	61.7,69.3,22,5,2.82,
	68.3,72.6,18,2,24.20,
	71.7,34.0,17,8,3.30,
	55.7,19.4,26,28,12.11,
	54.3,15.2,31,20,2.15,
	65.1,73.0,19,9,2.84,
	65.5,59.8,22,10,5.23,
	65.0,55.1,14,3,4.52,
	56.6,50.9,22,12,15.14,
	57.4,54.1,20,6,4.20,
	72.5,71.2,12,1,2.40,
	74.2,58.1,14,8,5.23,
	72.0,63.5,6,3,2.56,
	60.5,60.8,16,10,7.72,
	58.3,26.8,25,19,18.46,
	65.4,49.5,15,8,6.10,
	75.5,85.9,3,2,99.71,
	69.3,84.9,7,6,99.68,
	77.3,89.7,5,2,100.00,
	70.5,78.2,12,6,98.96,
	79.4,64.9,7,3,98.22,
	65.0,75.9,9,9,99.06,
	92.2,84.6,3,3,99.46,
	79.3,63.1,13,13,96.83,
	70.4,38.4,26,12,5.62,
	65.7,7.7,29,11,13.79,
	72.7,16.7,22,13,11.22,
	64.4,17.6,35,32,16.92,
	77.6,37.6,15,7,4.97,
	67.6,18.7,25,7,8.65,
	35.0,1.2,37,53,42.34,
	44.7,46.6,16,29,50.43,
	42.8,27.7,22,29,58.33),
	nrow = 47, ncol = 5, byrow = TRUE)

	y <- design[,1]
	x1 <- design[,2]
	x2 <- design[,3]
	x3 <- design[,4]
	x4 <- design[,5]

	model <- lm(y ~ x1 + x2 + x3 + x4)

	estimates <- matrix(c(91.05542390271397,6.94881329475087,
	-0.22064551045715,0.07360008972340,
	-0.26058239824328,0.27410957467466,
	-0.96161238456030,0.19454551679325,
	0.12441843147162,0.03726654773803),
	nrow = 5, ncol = 2, byrow = TRUE)

	expectedBeta <- estimates[,1]
	expectedErrors <- estimates[,2]

	expectedResiduals <- c(7.1044267859730512,1.6580347433531366,
	4.6944952770029644,8.4548022690166160,13.6547432343186212,
	-9.3586864458500774,7.5822446330520386,15.5568995563859289,
	0.8113090736598980,7.1186762732484308,7.4251378771228724,
	2.6761316873234109,0.8351584810309354,7.1769991119615177,
	-3.8746753206299553,-3.1337779476387251,-0.1412575244091504,
	1.1186809170469780,-6.3588097346816594,3.4039270429434074,
	2.3374058329820175,-7.9272368576900503,-7.8361010968497959,
	-11.2597369269357070,0.9445333697827101,6.6544245101380328,
	-0.9146136301118665,-4.3152449403848570,-4.3536932047009183,
	-3.8907885169304661,-6.3027643926302188,-7.8308982189289091,
	-3.1792280015332750,-6.7167298771158226,-4.8469946718041754,
	-10.6335664353633685,11.1031134362036958,6.0084032641811733,
	5.4326230830188482,-7.2375578629692230,2.1671550814448222,
	15.0147574652763112,4.8625103516321015,-7.1597256413907706,
	-0.4515205619767598,-10.2916870903837587,-15.7812984571900063)

	expectedStdError <- 7.73642194433223
	expectedRSquare <- 0.649789742860228
	expectedAdjRSquare <- 0.6164363850373927

	verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
	expectedStdError, expectedRSquare, expectedAdjRSquare, "Swiss Fertility")

	# model with no intercept
	model <- lm(y ~ 0 + x1 + x2 + x3 + x4)

	estimates <- matrix(c(0.52191832900513, 0.10470063765677,
	2.36588087917963, 0.41684100584290,
	-0.94770353802795, 0.43370143099691,
	0.30851985863609, 0.07694953606522),
	nrow = 4, ncol = 2, byrow = TRUE)

	expectedBeta <- estimates[,1]
	expectedErrors <- estimates[,2]

	expectedResiduals <- c(44.138759883538249, 27.720705122356215, 35.873200836126799,
	34.574619581211977, 26.600168342080213, 15.074636243026923, -12.704904871199814,
	1.497443824078134, 2.691972687079431, 5.582798774291231, -4.422986561283165,
	-9.198581600334345, 4.481765170730647, 2.273520207553216, -22.649827853221336,
	-17.747900013943308, 20.298314638496436, 6.861405135329779, -8.684712790954924,
	-10.298639278062371, -9.896618896845819, 4.568568616351242, -15.313570491727944,
	-13.762961360873966, 7.156100301980509, 16.722282219843990, 26.716200609071898,
	-1.991466398777079, -2.523342564719335, 9.776486693095093, -5.297535127628603,
	-16.639070567471094, -10.302057295211819, -23.549487860816846, 1.506624392156384,
	-17.939174438345930, 13.105792202765040, -1.943329906928462, -1.516005841666695,
	-0.759066561832886, 20.793137744128977, -2.485236153005426, 27.588238710486976,
	2.658333257106881, -15.998337823623046, -5.550742066720694, -14.219077806826615)

	expectedStdError <- 17.24710630547
	expectedRSquare <- 0.946350722085
	expectedAdjRSquare <- 0.9413600915813

	verifyRegression(model, expectedBeta, expectedResiduals, expectedErrors,
	expectedStdError, expectedRSquare, expectedAdjRSquare, "Swiss Fertility No Intercept")

	displayDashes(WIDTH)