Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE Estimation of Statistical Models in R Updated July 21, 2010 • Josh Mills Contents Introduction .............................................................................................................................................. 2 Importing Data ......................................................................................................................................... 3 Importing Data from Other Programs ......................................................................................................... 6 Ordinary Least Squares (OLS) Regression ................................................................................................... 9 Model Diagnostics ................................................................................................................................... 10 Two-Sample and K-Sample Tests (Univariate and Multivariate) ................................................................. 15 Heteroskedastic t Regression .................................................................................................................... 26 Fitting Data to Distributions .................................................................................................................... 27 Robust Regression ................................................................................................................................... 33 Poisson and Negative Binomial Models ..................................................................................................... 34 Zero-Inflated and Hurdle Poisson and Negative Binomial Models ................................................................ 37 Tobit Models, Sample Selection Models, and Truncated Models ................................................................. 39 Quantile Regression ................................................................................................................................. 43 Regression Scale Models .......................................................................................................................... 44 Partial Least Squares ............................................................................................................................... 46 Multivariate Adaptive Regression Splines .................................................................................................. 49 Regression Trees ..................................................................................................................................... 52 Nonparametric and Semiparametric Estimation ......................................................................................... 58 Nested Logit Models ................................................................................................................................ 61 Ordered Models ...................................................................................................................................... 65 Systems of Simultaneous Equations .......................................................................................................... 70 Seemingly Unrelated Regression ............................................................................................................... 74 Beta Regression ...................................................................................................................................... 78 Panel Data Models .................................................................................................................................. 79 Multilevel Models .................................................................................................................................... 85 Generalized Estimating Equations ............................................................................................................. 94 Time-Series Models ................................................................................................................................. 95 Using Google Maps in R .......................................................................................................................... 98 Using US Census Data with R .................................................................................................................. 99 Spatial Econometrics and Testing for Spatial Dependence ......................................................................... 103 Performing Analysis on Multiple Combinations of Files ............................................................................. 117 Estimation of Statistical Models Using R Page 1 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE Travel Demand Modeling in R ................................................................................................................ 123 Visualization of Networks and Graphs ...................................................................................................... 132 Large Data Sets: Virtual Memory and R .................................................................................................. 137 Miscellaneous Topics .............................................................................................................................. 140 Index of Packages Used .......................................................................................................................... 144 Introduction The R statistical package can best be thought of, in terms of user interface, as falling some- where in-between SAS and Limdep. Like SAS, R accomplishes different tasks using different procedures. Instead of the procedures being self-contained (like PROC REG in SAS), commands are issued one at a time, like Limdep. However, unlike SAS and Limdep, R can easily be extended using user-generated add-on packages, or libraries. Because R is open-source (and thus free), it is relatively easy to create add-on packages and quickly distribute them to a large, well-established community of users. Because of this, R can perform procedures and estimate statistical models that have not yet been implemented in purchaseable software such as SAS, Limdep, SPSS, and Stata. One primary example is spatial econometrics; while Stata contains capabilities for estimating certain types of spatial econometric models, the ma- jority of the most recent types of models (spatial panel models, spatial Probit, etc.) are only available in R. In short, R can give users access to cutting-edge methods and models easily. The available packages from the central R website, the Comprehensive R Archive Network (CRAN, http://cran.r-project.org), can accomplish s variety of tasks, such as statistical analysis (which is to be expected from a statistical software package), machine learning, importing and manipulating GIS data (R's GIS capabilities are far stronger than SAS), parallel compu- ting, qualitative research analysis, text mining, graph/network analysis, optimization, solving ordinary differential equations, image analysis, 3D graphics (through OpenGL), interactive graphics, and Bayesian inference. R can interface with a large number of other programs and platforms, such as Excel (through an add-in accessible through Excel), Google Maps, C++, WinBUGS, gretl, SAS, and even audio devices, relational databases (such as Oracle and SQL), the US Census database, LP_SOLVE, openNLP, PowerPoint, Word, LaTeX, the Web (through CURL), CPLEX, and Java. Despite this wide array of abilities, it is relatively easy to use a single data set with a wide variety of varieties, as will be demonstrated here. One important note - unlike SAS, all commands in R are case-sensitive. While package names and commands in the text of the handout are denoted using all caps, the exact capitali- zations should be used as shown in the syntax examples. Estimation of Statistical Models Using R Page 2 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE Importing Data Packages Used: base, foreign, utils, prettyR Use READ.TABLE or READ.CSV to import a data set (the CSV, or comma-separated values, format tends be more reliable than other formats such as Excel and Access, particularly in R, since there are no built-in functions for importing data from Excel or Access). In what may be surprising for a statistical software package, R can import data from URLs (in other words, a Web-based text file or data file can be imported through R simply by using a URL as the file name) and even directly from ZIP archives (which is accomplished using the unz command). READ.TABLE can import practically any delimited file (online or offline), but one must explicitly specify the record layout of the file, which can be tricky. In addition, importing files created by other statistical software packages, such as SAS, may lead to strangely specified variables. For example, SAS uses a period (".") to represent a missing value. When importing a file, R may misread this value and, instead of classifying a numerical variable as numerical (class NUMERIC), R would instead classify that variable as a factor (class FACTOR). Factors are useful for manipulating categorical data (for example, a variable in which a value of 1 indicates a rural minor arterial, 2 indicates a rural principal arterial, 3 indicates an urban minor arterial, etc.), but most data should be interpreted strictly as numeric values. To ensure this does not happen, use NA.STRINGS = "." (for data exported by SAS). Commands with similar syntax requirements that can import data files from Minitab, SAS, dBase, and Stata can be found in the FOREIGN package. For most commands, such as those involving importing data or estimating models, the results should be stored in separate variables. Unlike C++, one does not have to declare the "class" of a variable when it is created; R will automatically determine the appropriate class. Data imported using the READ.CSV command is stored in a data frame (class DA- TA.FRAME). Syntax trt<-read.csv("F:/lots_of_data/trtout.csv",header=F, na.strings = ".") If the first line of your data set has column names, set HEADER equal to T (for TRUE). In this example, the imported data will be stored in a data frame named TRT. Since the results of the command are being stored in the variable, one must refer to or access the variable in order to view the data or learn some characteristic about the data. To view the data, simply type the name of the variable and hit return. Syntax trt Output >trt V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 1 14 2 600 7 44 1 0 5 1 0 3 5 78 0 NA 2 13 2 400 7 44 0 0 1 0 1 0 1 85 1 NA Estimation of Statistical Models Using R Page 3 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE 3 11 2 540 7 44 1 0 2 1 0 0 2 81 0 NA 4 12 2 540 7 44 1 0 5 0 1 0 4 83 1 NA 5 NA 9 2 370 7 44 1 0 5 0 1 0 4 83 1 6 10 2 350 7 44 1 0 3 1 1 0 3 81 0 NA 7 10 2 350 7 44 1 0 3 0 1 0 1 80 0 NA 8 12 1 460 14 48 0 0 2 0 1 0 1 86 0 NA 9 NA 8 2 400 7 44 1 0 3 0 1 0 1 84 1 Note that because variable names were not provided, R has named each variable V1, V2, V3, and so on. Missing values are denoted by "NA." To verify that TRT is a data frame, use the CLASS command. Syntax class(trt) Output [1] "data.frame" It is easy to view only a subset of the rows or columns in the data frame. Use the following syntax to view the first three rows and the first four columns. Syntax trt[1:3,1:4] Output > trt[1:3,1:4] V1 V2 V3 V4 1 14 2 600 7 2 13 2 400 7 3 11 2 540 7 To refer to a particular variable, use the $ operator. In this example, a histogram for variable V9 will be displayed using the HIST command. Syntax hist(trt$V9) Output Estimation of Statistical Models Using R Page 4 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE R is particularly strong with graphics. The generated histogram can be easily pasted into a document or saved as an image file simply by right-clicking on the plot. Descriptive statistics can be viewed using the DESCRIBE and FREQ commands from the PRETTYR library. Syntax library(prettyR) describe(trt) freq(trt) Note that a library/package is loaded using the LIBRARY command. Output > describe(trt) Description of trt Numeric mean median var sd valid.n V1 13 13 7.4 2.72 111 V2 3.483 2 7.638 2.764 151 V3 364.9 385 7.116e+04 266.8 151 V4 135.9 10 5.643e+04 237.5 151 V5 37.93 43 414.5 20.36 151 V6 12.87 1 447 21.14 151 V7 0.3377 0 0.2918 0.5402 151 V8 1.987 2 1.786 1.337 151 V9 1.152 1 1.623 1.274 151 V10 0.6026 1 0.2411 0.491 151 V11 0.3311 0 0.3563 0.5969 151 V12 1.596 1 1.576 1.255 151 V13 60.45 81 1241 35.23 151 V14 22.43 1 1346 36.68 151 V15 0.525 1 0.2558 0.5057 40 frwy 0.06623 0 0.06225 0.2495 151 art 0.1589 0 0.1346 0.3668 151 frwytl 0.3642 0 2.006 1.416 151 rurtl 3.947 7 15.89 3.986 151 arttl 2.371 0 30.78 5.548 151 cage 25.55 5 1241 35.23 151 youngm 1 1 0 0 151 Estimation of Statistical Models Using R Page 5 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE speed 23.09 22 35.24 5.937 111 > freq(trt) Frequencies for V1 10 11 12 13 14 15 16 17 18 19 20 21 NA 32 6 12 14 21 9 5 4 3 2 1 2 40 % 21.2 4 7.9 9.3 13.9 6 3.3 2.6 2 1.3 0.7 1.3 26.5 %!NA 28.8 5.4 10.8 12.6 18.9 8.1 4.5 3.6 2.7 1.8 0.9 1.8 Frequencies for V2 1 2 3 6 7 8 9 NA 24 77 10 2 9 18 11 0 % 15.9 51 6.6 1.3 6 11.9 7.3 0 %!NA 15.9 51 6.6 1.3 6 11.9 7.3 … Variables can be added to a data frame by using the assignment operators (<--, -->, or =) and the $ operators. Some examples are listed below. Note the syntax of the IF-THEN statements. The "pipe" operator (|) is used for OR, and the double-ampersand (&&) is used for AND. Syntax trt$frwy<-ifelse(trt$V2==3,1,0) trt$art<-ifelse(trt$V2==1,1,0) trt$frwytl<-ifelse(trt$V2==3,trt$V4,0) trt$rurtl<-ifelse(trt$V2==2,trt$V4,0) trt$arttl<-ifelse(trt$V2==1,trt$V4,0) trt$youngm<-ifelse(trt$V9<3 && trt$V9 ==1,1,0) trt$speed<-(trt$V5/10)/(trt$V1/60) y<-trt$speed X<-trt[,1:ncol(trt)-1] The best way to learn about the syntax of a command is by using the HELP or ? command. Help files are provided in HTML format, making it easy to browse through the different commands (the index for the package being used can be found at the bottom of the help page. For example, the HELP command can be used to access the help file for the LM command in the STATS package, which is used to estimate models using ordinary least squares (OLS). Syntax help(lm) ?lm Importing Data from Other Programs Packages Used: foreign, xlsReadWrite, xlsx A number of add-on packages are available for importing data written by programs such as Microsoft Excel, SAS, SPSS, and Stata. R can also write to most (but not all) of these for- Estimation of Statistical Models Using R Page 6 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE mats. Excel files can be read and written using the XLSREADWRITE package. At the time of this writing, however, this package is not yet available for the 64-bit Windows version. The syntax is fairly straightforward, similar to using READ.CSV. Syntax ind<-read.csv("F:/Folder of Data/hoosierc.csv",header=T,na.strings=”.”) # Reading in Excel 2003 files using xlsReadWrite # Note that xlsReadWrite does not yet have a 64-bit version write.xls(ind,"F:/lots_of_data/hoosierexcel.xls",colNames=TRUE) indexcel<-read.xls("F:/lots_of_data/hoosierexcel.xls",colNames=TRUE) The XLSREADWRITE package can manipulate Excel 97-2003 files. The XLSX package can read and write files in the newer Excel 2007 (xlsx) format. The package, however, requires the installation of the Java Development Kit (http://java.sun.com/javase/downloads/index.jsp). This package can also directly manipulate Excel 2007 workbooks. The XLSX package uses Java-based commands. The documentation details only a subset of the operations that can be performed on Excel workbooks. Knowledge of Java is needed to understand the undocu- mented features. Data is read using the READ.XLSX command. This can be quite time-consuming for large files (> 100 columns). This is also true for the related WRITE.XLSX command. Syntax write.xlsx(ind,"F:/lots_of_data/hoosierexcel07.xlsx",sheetName="Indiana") # The write.xlsx command is quite time consuming indxlsx<-read.xlsx("F:/lots_of_data/hoosierc.xlsx",sheetName="hoosierc") # The read.xlsx command is also time-consuming for large files (> 1 MB) A suite of commands are available to perform commonly used tasks on Excel 2007 work- sheets, such as merging cells, changing print settings, and even adding pictures. A few of these are demonstrated below. Use help(loadWorkbook) to learn details about how these commands work. Syntax # The xlsx package also has Java-based commands to directly work with Excel spreadsheets # In this example, an xlsx file will be loaded and a picture will be added # The workbook will be set to print as landscape with 0.5 inch header and footer margins # The workbook will then be saved as a new file indwb <- loadWorkbook("F:/lots_of_data/hoosierc.xlsx") #.jmethods(indwb) # Java methods that can be used - must know Java! indsh <- getSheets(indwb) indsh1<-(indsh[['hoosierc']]) Estimation of Statistical Models Using R Page 7 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE indprint<-printSetup(indsh1, footerMargin=0.5,headerMargin=0.5, landscape=T,paperSize="LETTER_PAPERSIZE") addPicture("F:/lots_of_data/Indiana_state_seal.png", indsh1, star- tRow=10,startColumn=10) saveWorkbook(indwb,"F:/lots_of_data/hoosierp.xlsx") SAS, SPSS, Stata, and dBase (DBF) files can be read and written using the FOREIGN package. The DBF and SPSS formats are the most restrictive; both formats only allow column (variable) names of a maximum of eight characters. The syntax of this set of commands is straightforward. Syntax # Reading and writing SAS, SPSS, Stata, DBF formats using package "foreign" # DBF - dBase, FoxPro - Used with shapefiles # The DBF format is rather old - column names are limited to 8 characters # In this example, a CSV file will be read then exported out as a DBF # read.dbf will be used to import this file # R will give warning messages indicating truncation had occurred county<-read.dbf("F:/lots_of_data/co18_d00.dbf") write.dbf(ind,"F:/lots_of_data/hoosierdb.dbf") indf<-read.dbf("F:/lots_of_data/hoosierdb.dbf") # Reading in SPSS files inds<-read.spss("F:/lots_of_data/hoosierspss.sav",to.data.frame=TRUE) # Reading in and writing Stata (.dta) files # Version 11 can be obtained for a low price from Krannert, 7th floor write.dta(ind,"F:/lots_of_data/hoosierstata.dta",version=10) indstata<-read.dta("F:/lots_of_data/hoosierstata.dta") The WRITE.FOREIGN command can create program-specific importing code for SAS, Stata, and SPSS. In these examples, separate code files are generated for SAS (with extension .sas), Stata (as do-files, with extension .do), and for SPSS. Syntax # Using write.foreign to write data sets for SAS, Stata, and SPSS # This command will also write importing code for each respective program # SPSS export is picky - column names must be 8 characters or less foreign:::writeForeignSAS # View the underlying code write.foreign(ind,"F:/lots_of_data/hoosier1.csv", "F:/lots_of_data/hoosier11.sas", package=c("SAS"),dataname="hoosierR") write.foreign(ind,"F:/lots_of_data/hoosier2.csv", "F:/lots_of_data/hoosier22.do", Estimation of Statistical Models Using R Page 8 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE package=c("Stata")) write.foreign(ind,"F:/lots_of_data/hoosier3.csv", "F:/lots_of_data/hoosier33.txt", package=c("SPSS")) Ordinary Least Squares (OLS) Regression Packages Used: stats The LM command is used to estimate models by OLS. It is often convenient to store formu- las in a variable of class FORMULA to create more readable code. Syntax f1<-speed~1+frwy+art+cage+V6+V3 Instead of directly issuing the command, it is more beneficial to store the results of the LM regression in a variable and view the results using the SUMMARY command. Note that observations can be (temporarily) deleted using the NA.ACTION=NA.OMIT statement. Syntax f1<-speed~1+frwy+art+cage+V6+V3 ff1<-lm(f1,data=trt,na.action=na.omit) summary(ff1) Output > ff1<-lm(f1,data=trt,na.action=na.omit) > summary(ff1) Call: lm(formula = f1, data = trt, na.action = na.omit) Residuals: Min 1Q Median 3Q Max -10.3842 -2.6012 -0.4189 3.1523 14.5701 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 24.502845 1.927874 12.710 < 2e-16 *** frwy 9.371764 1.843557 5.084 1.62e-06 *** art 3.486268 1.292920 2.696 0.00817 ** cage 0.082596 0.110675 0.746 0.45716 V6 1.062831 0.989525 1.074 0.28525 V3 -0.008268 0.003288 -2.515 0.01343 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.885 on 105 degrees of freedom (40 observations deleted due to missingness) Multiple R-squared: 0.3537, Adjusted R-squared: 0.3229 F-statistic: 11.49 on 5 and 105 DF, p-value: 7.351e-09 Estimation of Statistical Models Using R Page 9 of 145 Institute of Transportation Engineers, Purdue Student Chapter engineering.purdue.edu/ITE The next section will describe a series of diagnostic tests that can be used to check the validity of the estimated model. Model Diagnostics Packages Used: car, lmtest, nortest, sandwich Compared to SAS and Limdep, a wider variety of statistical tests are available for testing the validity of estimated statistical models. A number of these tests can also be applied to generalized linear models (for Poisson and negative binomial models, for example). The CAR and LMTEST packages are used to conduct these tests. The SANDWICH package can be used to correct for heteroskedasticity. Four tests can be used to test for heteroskedasticity: the Score test, the Breusch-Pagan test, the Goldfeld-Quandt test, and the Harrison-McCable tests. For all of these tests, the null hypothesis is that the residuals of the estimated model exhibit uniform variance (no heteroskedasticity). Note the use of # to denote comments. Syntax ncv.test(ff1) # Score test for heteroskedasticity, car bptest(f1,data=trt) # Breusch-Pagan test for heteroskedasticity, lmtest gqtest(f1, data=trt) # Goldfeld-Quandt test for heteroskedasticity, lmtest > hmctest(f1, data=trt) > ncv.test(ff1) # Score test for heteroskedasticity, car Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 1.654428 Df = 1 p = 0.1983573 > bptest(f1,data=trt) studentized Breusch-Pagan test data: f1 BP = 5.3949, df = 5, p-value = 0.3696 > gqtest(f1, data=trt) Goldfeld-Quandt test data: f1 GQ = 1.0827, df1 = 50, df2 = 49, p-value = 0.3908 > hmctest(f1, data=trt) Harrison-McCabe test data: f1 HMC = 0.482, p-value = 0.403 Estimation of Statistical Models Using R Page 10 of 145

Estimation of Statistical Models in R

Similar Estimation of Statistical Models in R