{"id":314,"date":"2014-09-26T08:31:42","date_gmt":"2014-09-26T08:31:42","guid":{"rendered":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/?p=314"},"modified":"2025-02-26T13:21:38","modified_gmt":"2025-02-26T13:21:38","slug":"using-r-with-stata-part-ii","status":"publish","type":"post","link":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/09\/26\/using-r-with-stata-part-ii\/","title":{"rendered":"Using R with Stata: Part II"},"content":{"rendered":"

This is\u00a0the second\u00a0in a series of posting about conducting advanced statistical analyses in Stata by sending a job to R and then reading the results back into Stata. Our task for testing this process is to use the R package, DPpackage, to fit a Bayesian Dirichlet process mixture (DPM) model for smoothing a scatter plot. The data that we are using consist of 298 IgG measurements made on children of different ages.<\/p>\n

In my last posting<\/a> I described the R code needed to fit the model and this time\u00a0I will consider how that code can be submitted and run without leaving Stata.<\/p>\n

Let us first assume that we have a text file called DPC.R than contains the R script. Repeating what we had last time, the file will contain,<\/p>\n

setwd(‘E:\/IGG’)<\/span>
\n require(foreign)<\/span>
\n require(DPpackage)<\/span>
\n D <- read.dta(‘igg.dta’)<\/span>
\n zbar <- c(7,3)<\/span>
\n Sigma <- matrix(c(4,1,1,1),2,2)<\/span>
\n MyPrior <- list(alpha=5, k0=0.25, nu1=10, m1=zbar, psiinv1=7*solve(Sigma))<\/span>
\n MyDesign <- list(nburn=5000, nsave=5000, nskip=3, ndisplay=1000)<\/span>
\n xp \u2190 c(1,2,3,4,5)<\/span>
\n fitIgG <- DPcdensity(y=D$igg, x=D$age, xpred=xp,ngrid=100,<\/span>
\n compute.band=TRUE, type.band=”HPD”, prior=MyPrior, mcmc=MyDesign,<\/span>
\n state=NULL, status=TRUE)<\/span>
\n nClu <- as.data.frame(fitIgG$save.state$thetasave[,7])<\/span>
\n names(nClu) <- “nclus”<\/span>
\n write.dta(nClu,file=’nclusters.dta’)<\/span>
\n cDen <- as.data.frame(cbind(fitIgG$grid,t(rbind(fitIgG$densp.m,<\/span>
\n fitIgG$densp.l,fitIgG$densp.h))))<\/span>
\n names(cDen) <- c(“IgG”,”m1″,”m2″,”m3″,\u201dm4\u201d,\u201dm5\u201d,<\/span>
\n “l1″,”l2″,”l3″,”l4″,”l5″,”h1″,”h2″,”h3″,”h4″,”h5”)<\/span>
\n write.dta(cDen,file=’cDensity.dta’)<\/span><\/p>\n

Having created this text file, we can send it to R using Stata’s shell command. For instance, if we were to type the Stata command,
\n. shell notepad<\/span><\/p>\n

Stata would be suspended while the operating system opens the notepad text editor. Stata would become active again when notepad is closed.<\/p>\n

When R is installed from http:\/\/www.r-project.org\/<\/a> (R is free),\u00a0a bin directory\u00a0 is created that will include a file called R.exe. Within Stata, we can run\u00a0the R script from file DPC.R if we type,
\n. shell path<\/em>\/R.exe CMD BATCH DPC.R log.txt<\/span><\/p>\n

where the word path<\/em> is replaced by the full path to R.exe on your computer. This command will suspend Stata while R runs the script as a batch job (it might take a few minutes). R sends its output, including any error messages to the log file. We can subsequently inspect this log file within Stata using,
\n. type log.txt<\/span><\/p>\n

Now I run several different programs in this way and it is tedious to have to remember the path to each program so I store them in a file called executables.txt<\/strong>, which is saved in my PERSONAL directory.<\/p>\n

Here is a listing of the first 5 lines of my file executables.txt.<\/p>\n

PLINK,”C:\\Software\\PLINK\\plink.exe”<\/span>
\n FBAT,”C:\\Software\\fbat\\fbat202C_win\\fbat.exe”<\/span>
\n R,”C:\\Program Files\\R\\R-3.0.2\\bin\\x64\\R.exe”<\/span>
\n WINBUGS,”C:\\Software\\WinBUGS14\\WinBUGS14.exe”<\/span>
\n OPENBUGS,”C:\\Software\\OpenBUGS\\OpenBUGS.exe”<\/span><\/p>\n

The file contains the name of each program that I use, a comma, and then the path.<\/p>\n

Matched to this file I have an ado file called execute.ado<\/strong> that reads the required path from this file and runs the executable. The contents of execute.ado are<\/p>\n

program execute<\/span>
\n gettoken PG CD : 0<\/span>
\n local personal : sysdir PERSONAL<\/span>
\n tempname fh<\/span>
\n file open `fh’ using `personal’executables.txt ,read<\/span>
\n file read `fh’ inLine<\/span>
\n while !r(eof) {<\/span>
\n tokenize `”`inLine'”‘, parse(“,”)<\/span>
\n if upper(“`1′”) == “`PG'” local executable “`3′”<\/span>
\n file read `fh’ inLine<\/span>
\n }<\/span>
\n file close `fh’<\/span>
\n confirm file “`executable'”<\/span>
\n shell “`executable'” `CD’<\/span>
\n end<\/span><\/p>\n

So we could give Stata the command
\n. execute R CMD BATCH DPC.R log.txt<\/span>
\nand execute.ado would locate the appropriate path for R and \u2018shell\u2019 it together with the remainder of the command line.<\/p>\n

Our R script not only runs DPcdensity but it returns some of the output as Stata .dta files so after it has run we might plot a histogram of the number of cluster in each of the 5,000 stored iterations of the sampler.<\/p>\n

.\u00a0use nclusters.dta, clear<\/span>
\n. histogram nclus , start(0.5) width(1) freq title(Number of Clusters)<\/span><\/p>\n

\"r2fig1\"<\/a><\/p>\n

Alternatively we could plot the conditional distribution at age 5 years using,<\/p>\n

.\u00a0use cDensity.dta, clear<\/span>
\n. twoway (line m5 igg , lpat(solid) lcol(red)) \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (line l5 igg , lpat(dash) lcol(black)) \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (line h5 igg , lpat(dash) lcol(black)), \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 leg(off) title(Conditional distribution of IgG) ytitle(f(IgG|Age=5))<\/span><\/p>\n

\"r2fig3\"<\/a><\/p>\n

This method can be adapted for any R analysis but it has one severe limitation. To use it, you must be able to program in R. This is fine for my own work, but is not very helpful if I am collaborating with someone who is not an experienced R programmer.<\/p>\n

Next time I will describe how I make the facilities of R available within Stata for people who do not know how to program R. The approach involves writing a Stata ado file that creates the R script as well as running it. In this way a Stata user only needs to install R on their computer, they do not need to know anything of its programming language.<\/p>\n","protected":false},"excerpt":{"rendered":"

This is\u00a0the second\u00a0in a series of posting about conducting advanced statistical analyses in Stata by sending a job to R and then reading the results back into Stata. Our task for testing this process is to use the R package, DPpackage, to fit a Bayesian Dirichlet process mixture (DPM) model for smoothing a scatter plot. […]<\/p>\n","protected":false},"author":134,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[37],"class_list":["post-314","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-r-stata-shell-executables-txt"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/314","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/users\/134"}],"replies":[{"embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/comments?post=314"}],"version-history":[{"count":8,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/314\/revisions"}],"predecessor-version":[{"id":338,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/314\/revisions\/338"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=314"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=314"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=314"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}