A Population Regression

Here's a video on some of the theory behind simple linear regression.



There's no R involved with this video, but the video provides some theory behind what it is that R's lm() command estimates.

Video Tutorial on Robust Standard Errors

Update: I have included a modified version of this summaryR() command as part of my package tonymisc, which extends mtable() to report robust standard errors. The tonymisc package is available on CRAN through the install.packages() command.

If you have the right R commands at your disposal, it is simple to correct for heteroskedasticity using the robust correction that is commonly-used among economists. I recorded a video tutorial to describe the simplest (and most flexible) way I know to get R to compute robust standard errors.

The key is to use a "summary-style" command that has an option to correct for heteroskedasticity. The command I like to use is called summaryR(). Here is the script file with the summaryR() command.



I found this function on an R-help discussion board where several people were answering someone's question about extending the summary.lm() command.

I deserve none of the credit for writing this (credit goes to John Fox), but I consider it my duty to point out how nice this function is. I demonstrate how to use the function in this video



Here are is the script file I used in the video:



Here's a link to the data set.

Video Tutorial on IV Regression

Update: [1 May 2011] I am working on a better augmentation of the current IV regression functions (specifically ivreg() in AER) in R. I will post a link here to my new method/function for IV regression when I finish debugging the code.

Update 2: [15 May 2011] The function I wrote here does not work with factors due to the way in which I constructed fitted values for the standard error correction in the second stage.

I recorded a new video tutorial on how to conduct an IV regression. This time I used my new ivregress() command, which has much better syntax than my old one. There are several reasons I needed to do this:

1. Relative to my previous post on IV regression, I have added the ability to conduct overidentifying tests. Now, my ivregress() command has all of the basic functionality of Stata's 2sls option to ivregress.
2. Relative to my previous video tutorial on IV regression, this video uses my newer command, which as much better syntax. As such, I will use this video tutorial to replace the old one.
   
   

Here is the R script file I used to define my ivregress() command and its companion summary command sum.iv().



To get ivregress() and sum.iv() to work, (1) copy the text in the above text box into an R script, (2) run all (or just download and run the R script) and (3) save into your default workspace if you want to use the commands into the future. Then, you're good to go. See my video demonstrating how to use the command here:



Here is the code I used to produce the output in the video.



Finally, here is a link to the synthetic "market share" data I used in the video.

I hope you find this command helpful.

An easier to use IV regression command in R

Update: I have added some functionality to my ivregress() command. Check out my newer post here.

After I posted my last video tutorial on how to use my IV regression function, I received a comment asking why I didn't write the command a different way to make the syntax easier to read.

The answer is that I didn't know how to write an easier to use function a year ago (when I wrote the ivreg() function). After some digging, I figured out how to work with "formula objects" in R and the result is an easier to use IV regression function (called ivregress()).

How to "install" ivregress()

Here's the code you need to run to define ivregress() and its companion summary command sum.iv(). This will provide instrumental variables regression estimates if you have one endogenous regressor with one or more instruments.



You only need to run the above code once to define the function object ivregress() for all future uses (as long as you don't write over it or clear your workspace).

How to use ivregress()

To create an "IV object" simply use the ivregress command as follows (the command relies on the car library):

library(car)
myivobject = ivregress(Y~X1+X2+...+Xk, Xj ~ Z1+Z2+...+Zl, dataframe)

where X1, X2,...Xk are second stage regressors, Xj is an endogenous regressor for which you would like to instrument, and Z1, Z2,...,Zl are instruments.

Then use sum.iv() from the previous posts to produce summary output. Sometime soon, I will post a video tutorial showing how to use ivregress() to perform IV regression.

Video Tutorial on Instrumental Variables in R

Update: I have replaced this video tutorial with a video tutorial on a newer, easier to use IV regression command. Check out that command here.

In this video, I show how to use my instrumental variables function in R, ivreg(), along with its companion summary command, sum.iv().



Here is a link to the synthetic "market share" data I used in the video.

IV Regression

Update: Since this post, I discovered a bunch of better ways to conduct 2SLS regression in R. These methods do not allow the user to look under the hood to see what is going on (so this post may be instructive anyway). If you're interested, ivreg() in AER, tsls() in sem and systemfit() in systemfit are good options. My package tonymisc also has a command, iv() that works when some or all of these fail. I'm working on making that function better (5/18/2011)

Here is my code from a previous post that performs IV regression. This may be easier to copy into an R script. I will post a video tutorial using this code shortly.

How Does R Deal with Categorical Explanatory Variables?

Here's a video tutorial where I demonstrate an answer to this question.



This video is designed to instill a baseline level of practical knowledge. There is more to how R treats factors in regression models. An interested reader should Google "R contrasts" for more information.

Here is the code I used in the video:



Here is a link to the data set.

Running a Regression in R

I created another video tutorial on R. This time, I discuss R's lm() command and how to use it for a variety of standard applications.



Here is the code that goes with the video:

Enjoy!

Two-stage least squares

Update: Since this post, I discovered a bunch of better ways to conduct 2SLS regression in R. These methods do not allow the user to look under the hood to see what is going on (so this post may be instructive anyway). If you're interested, ivreg() in AER, tsls() in sem and systemfit() in systemfit are good options. My package tonymisc also has a command, iv() that works when some or all of these fail. I'm working on making that function better (5/18/2011)

Here is some code I wrote that allows you to easily perform two-stage least squares regression using R.

ivreg = function(Y0,X0,Z10=NULL,Z20,data.df){   
  Y  = data.df[[Y0]]
  X  = data.df[[X0]]
  Z1 = NULL
  if(length(Z10)>0){
    Z1 = data.df[[Z10[1]]]
  }
  Z2 = data.df[[Z20[1]]]

  if(length(Z10)>1){
     for(i in 2:length(Z10)){
        Z1=cbind(Z1,data.df[[Z10[i]]])
     }
  }

  if(length(Z20)>1){
     for(i in 2:length(Z20)){
        Z2 =cbind(Z2,data.df[[Z20[i]]])
     }
  }

## First Stage
  Z       = cbind(1,Z1,Z2)        

  zPz     = t(Z)%*%Z
  zPx     = t(Z)%*%X
  zPz.inv = solve(zPz)
  beta.f  = zPz.inv%*%zPx

  xhat    = Z%*%beta.f
  SSE.f = sum((X - xhat)^2)
  df.d    = length(data.df[,1])-length(Z[1,])
  MSE.f   = SSE.f/df.d
  RMSE.f  = sqrt(MSE.f)
  se.f    = RMSE.f*sqrt(diag(zPz.inv))
  t.f     = beta.f/se.f
  pval.f  = 2*pnorm(abs(t.f),lower.tail=F)

## Constructing F-test
  Z.d       = cbind(rep(1,
                       length(data.df[,1])),Z1)
  zdPzd     = t(Z.d)%*%Z.d
  zdPx      = t(Z.d)%*%X
  zdPzd.inv = solve(zdPzd)
  bhat.d    = zdPzd.inv%*%zdPx
  xhat.d    = Z.d%*%bhat.d

  SSE.d      = sum((X-xhat.d)^2)
  SSEx.z     = SSE.d -SSE.f
  MSEx       = SSEx.z/length(Z20)
  F.stat     = MSEx/MSE.f
  df.n       = length(Z[1,]) -
                    length(Z.d[1,])
  pval.Ftest = pf(F.stat, df.n, df.d ,
                lower.tail=FALSE)

  ftest  = round(c(F.stat,
                   pval.Ftest,df.n), 4)
  names(ftest) = c("F-stat", "P-value",
                  "# of instruments")
  rn.f = c("constant",Z10,Z20)

  first   = data.frame(cbind(beta.f,se.f,
                      t.f, pval.f),row.names=rn.f)
  names(first) = c("Estimate", "Std. Error",
                  "T-stat", "P-value")

## Second Stage
  X.mat   = cbind(1,Z1,xhat)
  X.mat2  = cbind(1,Z1,X)

  xPx      = t(X.mat)%*%X.mat
  xPx.inv  = solve(xPx)
  xPy      = t(X.mat)%*%Y
  beta2sls = xPx.inv%*%xPy

  resid    = Y-X.mat2%*%beta2sls
  MSE      = mean(resid^2)
  RMSE     = sqrt(MSE)
  se2sls   = RMSE*sqrt(diag(xPx.inv))
  tstat    = beta2sls/se2sls
  pval     = 2*pnorm(abs(tstat),lower.tail=F)

  rn = c("constant",Z10,X0)

  second   = data.frame(cbind(beta2sls,se2sls,
                      tstat, pval),row.names=rn)
  names(second) = c("Estimate", "Std. Error",
                    "T-stat", "P-value")
  result = list(first,second,ftest)
  names(result) = c("first", "second", "ftest")

  print("IV object successfully created.
           Use sum.iv() on object")
  print("to learn about your 2SLS
         Regression")
  return(invisible(result))

}

The above function is ivreg(), which takes five arguments: (1) Y0 is the name of the dependent variable, (2) X0 is the name of the endogenous regressor, (3) Z10 is the set of names of exogenous regressors [second stage variables assumed to be exogenous], and (4) Z20 is the set of instruments.

The function will not print any output if run alone. The appropriate syntax is to save the function output into an "iv" object, which you can summarize using the following function:

sum.iv = function(reg.iv, first=FALSE,
        second=TRUE, ftest=FALSE) {
x= rep(0,2)
if(first==TRUE) x[1] = 1
if(second==TRUE) x[2]= 2
if(ftest==TRUE) x[3]= 3

print(reg.iv[x])
}

Here is an example of how to use the functions to perform a two-stage least squares regression. Imagine that you have a data on prices and quantities of cars sold, and you have some data on attributes that directly affect quantity sold (falls into category Z10) and attributes that only affect quantity sold through the price (falls into category Z20).

You could estimate the model, and summarize its output using the following two commands:

demand.iv = ivreg("quantity", "price", c("horsepower", "mpg"), c("regulations", "cost", "competition"))

sum.iv(demand.iv, first=FALSE, second=TRUE, ftest=FALSE)

Note: The quotes in the ivreg command are necessary to get the function to work properly. This is an unusual artifact of how I wrote the function.

What is nice is that the sum.iv() command allows the user to specify what output he/she wants to see from the two-stage least squares regression. For the sum.iv() command, you could just leave off the logical arguments to get the same output (because these are the defaults), but you can choose to show the first-stage estimates and the F-test on the explanatory power of the instruments by specifying your own options.