![\"nb_cov_test\"](\"http:\/\/blogs.oregonstate.edu\/geo599spatialstatistics\/files\/2014\/05\/nb_cov_test.png\")
<\/a><\/p>\nTo make sure that we only compare variables of interest, I preface the ggpairs() call with a few lines of code that build a data.frame containing only the variables we want to compare.\u00a0 The line “#var_tab$Acreage <- as.integer(var_tab$Acreage)” begins with a #, which is how you specify comments in R.\u00a0 That means this line is “commented out”.\u00a0 We will be running a negative binomial regression, and R assumes you will provide count data as the dependent variable.\u00a0 The regression still operates with float data, but R will issue warnings, which you can safely ignore.\u00a0 If you do not want to see warnings, simply round off all the decimals by uncommenting this line and running it.<\/p>\n
Note that I have commented out another line in the script: #ct <- qt(.95, df=(nrow(var_tab) – length(independentVars))).\u00a0 I didn’t end up using this because the package I used to report coefficient results did not allow me to specify different confidence intervals using critical values, but this is too bad.\u00a0 I do not have enough farms to justify a .95 probability and should be using a t distribution to determine confidence intervals.<\/p>\n
You may be wondering: why negative binomial regression?\u00a0 Well, “number of acres” is essentially count data.\u00a0 I have no reason for assuming it to be normally distributed around any particular mean acreage.\u00a0 In fact, I expect smaller farms to be more common than larger, so it is not unreasonable to hypothesize that the distribution of farms would match a poisson or negative binomial distribution.\u00a0 Turns out, the poisson regression was over-dispersed, so the negative binomial is the appropriate distribution in this case.\u00a0 It seems to me that the poisson is always over-dispersed, and I am getting into the habit of skipping straight to negative binomial.\u00a0 For an interesting analysis of the suitability of different distributions for analyzing ecological data, I recommend a paper by Walter Stroup called “Rethinking the Analysis of Non-Normal Data in Plant and Soil Science”.<\/p>\n
Go ahead and run the regression and examine the results using the command “summary()”.\u00a0 Note that the spatial component of my analysis “miles”, which represents the distance of the farm from the Corvallis Farmer’s Market in miles, is not significant.\u00a0 This is a shame because I hypothesized that farms would become smaller the closer they were to the market.\u00a0 In particular, I believe that there is a new business model for farms emerging based on small acreages producing high quality food intended for direct-to-consumer local markets.\u00a0 These results show that older farms tend to have larger acreage, and the higher proportion of their production that goes to local markets, the smaller the acreage tends to be.<\/p>\n
The package texreg is for formatting model results into the LaTeX language, which many academics use for publishing because it specializes in formatting model formulas and symbols.\u00a0 The function “plotreg()” makes an attractive graphic of the coefficient results.\u00a0 The next few lines of code simply print a Q-Q plot and a plot of the residual vs. fitted values, so we can test the model for normality and independence.\u00a0 The errors are reasonably normally distributed.\u00a0 You can see in the residuals vs. fitted plot that most of the farms are very small in size, but the errors seem independent of the farm size.\u00a0 Hard to tell with so few data points!<\/p>\n Next we append the residuals and how many standard errors each residual deviates from the mean to the shapefile data.\u00a0 This is the same syntax used to add data to a data.frame.\u00a0 The function writeSpatialShape() resides in the maptools package and does exactly what it sounds like, depositing the new shapefile in the working directory.\u00a0 Then, since we cannot see the output from R when calling the script remotely in ArcGIS, the code prints the results in a .dbf file.\u00a0 The function write.dbf() is in the package foreign.<\/p>\n Information criteria are all the rage now, so the final part of the code makes a .dbf table that includes every possible combination of independent variables included in the original model fit against the dependent variable, with the resulting AIC and BIC.\u00a0 Unless you have an abiding fascination with for and if loops, you can safely ignore the algorithm I use to permutate through every unique combination of independent variables.\u00a0 It stores all the combinations as vectors of numbers in a list.\u00a0 The code after it connects the variable names to the numbers and then pastes them together into formulas.\u00a0 This part of the code is only interesting if you work with list type data and want to know how to extract and manipulate data contained in lists.\u00a0 Hint: use the command “unlist()”.<\/p>\n Use write.dbf again to produce a table you can load in ArcGIS, and feel free to run the line containing just the name of the table “diag_tab” to view the diagnostics table in R. This concludes a poor tutorial in R.<\/p>\n If you are interested in how ArcGIS and R can talk to each other, here are a few tips to get you started.<\/p>\n I found this video by Mark Janikas<\/a> describing a python script that calls R to perform a logit regression and return the results to a shapefile in ArcGIS.\u00a0 The video is not super helpful in explaining how this process works, its intended audience is a conference full of pro programmers who are hip to hacking lingo, which is not me.\u00a0 At the end, however, he refers the interested to this github site<\/a> that hosts the python and r scripts.<\/p>\n One important note: certain ArcGIS directories have to be added to the search path for Python.\u00a0 You can manually add these directories to the search path by modifying the path variable in the environmental variables according to the instructions in the folder demo_geo599\/R_Tools\/Doc.\u00a0 However, I did this and it did not work, so I used a workaround.\u00a0 I manually add the paths using the code in the python file “pathlib.py” which I include in the demo folder.\u00a0 This file and the file “arcpyWithR.py” need to be in the python search path somewhere for my example file to work.\u00a0 And yes it will work on your home machine, just not in the lab.\u00a0 I have added a python file and R file titled “GLMwithR” to the folder R_Tools\/Scripts.<\/p>\n If you open the python file, you will note the first line reads “import pathlib”, which adds the necessary directories to the search path.\u00a0 Scrolling down, you will see the function “ARCPY.GetParameterAsText()” called six times.\u00a0 Each of these commands receives user input through the ArcGIS GUI interface.\u00a0 You can add the toolbox inside the R_Tools folder to your ArcGIS toolbox by right-clicking on the directory in the toolbox, selecting “Add Toolbox…” and navigating to the “R Tools” toolbox in the R_Tools folder.\u00a0 I added my script to this toolbox by right-clicking on the box and selected “Add -> Script…”.\u00a0 I would refer you to the help file on creating script tools<\/a>.\u00a0 The hardest part is actually configuring the parameters.\u00a0 If you right-click on the “Negative Binomial GLM” tool and select “Properties”, you can view the parameters by selecting the “Parameter” tab.<\/p>\n Note that there are six parameters.\u00a0 These parameters need to be specified in the same order as they appear on the python script, and each parameter needs to be configured properly.\u00a0 In the Property field, you have to select the proper values for Type, MuliValue, Obtained from, etc., or the script will fail to work.\u00a0 It’s pretty annoying.<\/p>\n I don’t understand all the python code myself, but the middle portion passes all the inputs to R in the list called “args”.\u00a0 Feel free to open the r script “GLMwithR.r” in the folder R_Tools\/Scripts and note the differences between it and the tutorial file we just worked through.<\/p>\n Python uses the functions “ARCPY.gp.GetParameterInfo()” and “ARCPY.SetParameterAsText()” to read the new shapefile and the two .dbf tables respectively that come from the R script.\u00a0 The whole process of sending info in between Python and R is still mysterious to me, so I can’t really shed any light on how or why the python script works.<\/p>\n That concludes my trail of breadcrumbs.\u00a0 Good luck!<\/p>\n <\/p>\n <\/p>\n","protected":false},"excerpt":{"rendered":" I would like to show you how to use a script tool within the ArcGIS toolbox to call R remotely, have it perform some function on your data, and return the results into a new shapefile projected onto your workspace in ArcMap.\u00a0 I won’t be doing that. The scripting tool works.\u00a0 I have used it… Continue reading <\/a><\/p>\n<\/a> <\/a><\/p>\n<\/a><\/p>\n<\/a><\/p>\n