## To run this demo, you will need to have downloaded R from
## http://cran.r-project.org/ and started an R session.  Type in the
## commands below at the command prompt (without the prompt itself and 
## without the comments that follow a # symbol), and type <Enter> after 
## each command. 

> ?cars             # Gives information about the data.  Press 'q' to exit this
  		      screen.

> cars              # Displays the data

> plot(cars,xlim=c(0,25))        # Keep this new window open

> X <- cars[,1]     # Set X to be the first column of the cars data
> X                 # Note X contains the car speeds
> n <- length(X)

> Y <- cars[,2]  
> Y

> beta.hat <- solve(t(X) %*% X) %*% t(X) %*% Y   # 'solve' computes a matrix 
# inverse; t(X) is the transpose of X; %*% denotes matrix multiplication

> beta.hat

> ?abline                   # To find out about the abline function
> abline(a=0,b=beta.hat)    # abline adds a line with intercept a and slope b 

## Now let's try using R to fit a linear model

> Mod1 <- lm(Y~X)           # Fits the model Y_i = alpha + beta X_i + epsilon_i
                            # The intercept is automatically included 
> names(Mod1)               # The output of the lm function is an object with 
                            # many names.  

> Mod1$coef                 # Extracts the estimated coefficients
> Mod1$coef[1]              # The first component
> Mod1$coef[[1]]            # Just the numerical value

> abline(Mod1)              # Notice how the abline function can be applied to 
                            # a linear model object        

## The fitted lines are different because our original model did not
## contain an intercept term.  Let's see if we can reproduce the
## calculation performed by the lm function.

> ?cbind
> X0 <- cbind(rep(1,n),X)  # Creates a new design matrix
> X0 

> beta.hat0 <- solve(t(X0) %*% X0) %*% t(X0) %*% Y 
                            # What is being computed here?
> beta.hat0             # Check this agrees with the coefficients of Mod1

## Theory suggests we should include a speed^2 term in our model

> Xstar <- cbind(X,X^2)  
> Xstar 

> Mod2 <- lm(Y~Xstar - 1) # Fits Y_i=alpha + beta X_i +gamma X_i^2 + epsilon_i 
                           # Notice that this is still a linear model 
> Mod2$coef[[1]]

> x <- seq(0,30,length=151) 
> x
> plot(cars,xlim=c(0,25))                # Replot the data without the lines
> lines(x,Mod2$coef[[1]]*x + Mod2$coef[[2]]*x^2)

## Exercise: What commands are needed to reproduce the calculation
## performed by the lm function to compute the estimated coefficient
## vectors in Mod2?