Parameter Estimation (Reading: Faraway (2005, 1st edition), section 2.8)

Here is a data set concerning the number of species of tortoise on the various Galapagos Islands. There are 30 cases (Islands) and 7 variables in the data set. We start by reading the data into R.

> gala <- read.table("gala.data") # read the data into R
> gala # take a look

The variables are

In this dataset,

• Some missing values are filled in for simplicity.

• Q: Are these variable quantitative or qualitative? If quantitative, continuous or discrete?

Now fit a model. Let us start with a simple model, in which

• the variable "Species" is treated as the response (Q: Is it reasonable to regard Species as a continuous variable?) and

• the rest variables except "Endemics" as the predictors.

• The model can be expressed as

The command to fit a linear model in R is:

> gfit <- lm(Species~Area+Elevation+Nearest+Scruz+Adjacent, data=gala) # "lm" is the command in R to fit a linear model. Use "help(lm)" to understand the details of the command. The option "data=" specifies the dataset that will be used. Also, use "help(formula)" to learn how to express a model in R.
> summary(gfit) # all the usual regression stuff

 Notice that

• There are many quantities (statistics) and many information in the output.

• Q: What quantities have caught your attention?

• For different datasets, you might find their most important/useful information existing in different statistics.

• Q: For this dataset, which statistics/quantities are particularly important?

The default of "summay(lm.object)" does not print the correlation matrix.

> summary(gfit, cor=T) # with the option "cor=T", the output will contain the correlation matrix of the estimated parameters

> options(digits=3) # the option "digits" controls the number of digits to print

> summary(gfit, cor=T)$cor # take a look of the correlation matrix of β's estimators

Notice that

• the correlation matrix is a symmetric matrix

• its diagonals are 1's

• Q: How to read these numbers?

  • For example, cor(β₁-hat,β₂-hat)=-0.8014 shows that the estimators of Area and Elevation effects are highly negatively correlated. They have strong collinearity.

  • On the other hand, cor(β₂-hat,β₄-hat)=0.099 tells us the estimators of Elevation and Scruz effects have much weaker correlation.

> summary(gfit)$cov * (summary(gfit)$sigma^2) # the covariance matrix of the β's estimators, the command "summary(gfit)$cov" returns (X^TX)^-1

Notice that in the covariance matrix

• The values in the diagonal are the estimated variances of β-hat. Compare it with the Std. Error in the lm output. Q: What is the relationship between them?

• Q: How are the covariance matrix related to the correlation matrix? In other words, how to obtain the correlation matrix from the covariance matrix given above?

Let us now examine the correlation between the columns of X.

> cor(model.matrix(gfit)) # the command "model.matrix" returns the model matrix X, which in this case is exactly the matrix formed by the data of the five predictor variables themselves plus the constant vector. This output shows the correlation between the columns of X.

Notice that

• this correlation matrix is related to X^TX, rather than the (X^TX)^-1 in the correlation matrix of β-hat.

• Check how different it is from the correlation matrix of β-hat. Q: Can you identify some relationship between the two correlation matrices?

• Q: Why do we get the NA's in the output?

> cov(model.matrix(gfit)) # the covariances between the five variables.

Check how different it is from the covariance matrix of the estimated parameters.

For practice purpose, let us also construct these quantities directly --- we define Y, X, (X^TX)^-1, %*% denotes matrix multiplication, "t()" makes the transpose, and "solve()" computes the inverse.
> y <- gala$Species # the response
> rep(1,30) # make a vector of ones

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

> x <- cbind(rep(1,30), gala[,-c(1,2)]) # bind a column of ones to all except the first two columns of the data from gala
> x <- as.matrix(x) # force x from being a data frame to a matrix x
> xtxi <- solve(t(x) %*% x) # compute (X^TX)^-1
> xtxi # it is the same as the output obtained from the command "summary(gfit)$cov"

β-hat, y-hat, ε-hat may all be computed directly or extracted from the linear model object gfit. We do it both ways here --- check whether the results match.

> beta <- xtxi %*% t(x) %*% y # the β-hat
> beta 

However,

• it is a bad way to compute β-hat

• it is inefficient and can be very inaccurate when the predictors are strongly correlated

• An alternative is to directly solve the normal equation X^TXb=X^TY.

  > solve(t(x)%*%x, t(x)%*%y) # solve the normal equation

                 [,1]

  rep(1, 30)  7.06822

  Area       -0.02394

  Elevation   0.31946

  Nearest     0.00914

  Scruz      -0.24052

  Adjacent   -0.07480

We can extract the regression quantities we need from the lm object.

> names(gfit) # the "names()" command is the way to see the components of an R object

> gfits <- summary(gfit)
> names(gfits) # more regression quantities in summary of lm object.

> deviance(gfit) # the "deviance" command returns RSS for a regression fitted model.

[1] 89231

Compare the quantity and the Residual standard error in the lm output. Q: how to obtain the residual standard error from the RSS? Notice that σ-hat=sqrt(RSS/(n-p)).

> gfit$coef # get β-hat from the fitted model gfit

> yhat <- x %*% beta # y-hat=X(β-hat)
> cbind(yhat, gfit$fitted) # compare if they are identical

> res <- y - yhat # the residuals
> cbind(res, gfit$resid) # compare if they are identical

Notice that

• the residuals and the fitted values are orthogonal.

> sum(res*yhat)  # their inner product = 0

[1] 2.73e-11

• actually, the residual vector is orthogonal to any columns of X (therefore, orthogonal to any linear combinations of the columns of X

  > t(x) %*% res # X^T(ε-hat)

                  [,1]

  rep(1, 30)  4.76e-12

  Area        1.54e-09

  Elevation   2.91e-11

  Nearest    -6.96e-13

  Scruz       6.98e-11

  Adjacent   -5.07e-10

We now compute the residual sum of squares (RSS), σ-hat, standard error of β-hat, R², and the correlation matrix for β-hat. Verify that the results agree with the output from before.

> sum(res*res) # the RSS

[1] 89231

> gfit$df # the degrees of freedom, it = 30-6, why?

[1] 24

> sighat <- sqrt(sum(res*res)/gfit$df) # calculate σ-hat
> sighat # take a look

[1] 61

> diag(xtxi) # the "diag" command extract the diagonal of a matrix. The output is the diagonal entries of (X^TX)^-1

> sqrt(diag(xtxi))*sighat  # Q: what's this?

> summary(gfit)$coef[,2]  # compare with the previous outputs

> mean(y) # calculate the mean of y

[1] 85.2

> sum((y-mean(y))^2) # the total sum of squares

[1] 381081

> 1 - sum(res*res) / sum((y-mean(y))^2) # the R²

[1] 0.766

> var(yhat)/var(y) # Q: why does it equal R²?

> cor(y,yhat)^2 # Q: why does it equal R²?

Let's calculate the correlation matrix of β-hat.

> z <- sqrt(diag(xtxi))
> xtxi / z %o% z # the correlation matrix, where %o% means outerproduct zz^T

Now we make a prediction for the number of species on an island with predictor values: Area=0.08, Elevation=93, Nearest=6.0, Scruz=12.0, and Adjacent=0.34.

> xp <- c(1, 0.08, 93, 6.0, 12.0, 0.34) # don't forget the initial one
> sum(xp*beta) # the predicted number of species