my last posting<\/a> I described how we can use Stata\u2019s \u2013shell- command to submit an R job and this time we will consider how\u00a0the R script\u00a0can be\u00a0created by\u00a0a Stata ado file so that a Stata user need know nothing of R when they fit the model.<\/p>\nPreviously, we worked with a text file called DPC.R than contains the R script. Repeating what we had last time, the file contains,<\/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 <- 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 head(cDen)<\/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>\nLast time we prepared this script, submitted it from within Stata and read the results into Stata for further processing. This approach is fine for those of us who know R but let’s now suppose that we want to prepare a program for those Stata users who know nothing of R. We will need to prepare a program that writes the script for them as well as executing it.<\/p>\n
There are lots of ways of going this but I will describe my preferred method.<\/p>\n
Within my PERSONAL folder I have a subfolder called RTEMPLATES and here I store the templates for the R scripts that run the jobs that I want to make available. Here is the template corresponding to DPC.R.<\/p>\n
setwd(‘@’)<\/span>
\nif( !require(‘foreign’) ) {
\ninstall.packages(‘foreign’,repos=’http:\/\/cran.ma.imperial.ac.uk\/’,dep=TRUE)
\nif( !require(‘foreign’) ) stop() }
\nif( !require(‘DPpackage’) ) {
\ninstall.packages(‘DPpackage’,repos=’http:\/\/cran.ma.imperial.ac.uk\/’,dep=TRUE)
\nif( !require(‘DPpackage’) ) stop() }<\/span>
\nD <- read.dta(‘@’)<\/span>
\n zbar <- c(@)<\/span>
\n Sigma <- matrix(c(@),2,2)<\/span>
\n MyPrior <- list(alpha=@, k0=@, nu1=@, m1=zbar, psiinv1=@*solve(Sigma))<\/span>
\n MyDesign <- list(nburn=@, nsave=@, nskip=@, ndisplay=1000)<\/span>
\n xp <- c(@)<\/span>
\n fitIgG <- DPcdensity(y=D\u00a3@, x=D\u00a3@, xpred=xp, ngrid=@,<\/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\u00a3save.state\u00a3thetasave[,7])<\/span>
\n names(nClu) <- “nclus”<\/span>
\n write.dta(nClu,file=’nclusters.dta’)<\/span>
\n cDen <- as.data.frame(cbind(fitIgG\u00a3grid,t(rbind(fitIgG\u00a3densp.m,<\/span>
\n fitIgG\u00a3densp.l,fitIgG\u00a3densp.h))))<\/span>
\n names(cDen) <- c(@)<\/span>
\n write.dta(cDen,file=’cDensity.dta’)<\/span><\/p>\nA template is almost the same as the script itself but with three minor changes; the symbol @ replaces items that are problem dependent, such as the values of the priors, these will be replaced\u00a0before the program is run;\u00a0the symbol \u00a3 replaces $ and there is a little extra code associated with loading packages. The \u00a3’s stop Stata from treating lines containing code, such as D$age, as if they contain a Stata global and the\u00a0extra code associated with loading packages\u00a0is intended to protect the Stata user from the need to install the package.<\/p>\n
In R, packages are installed over the internet and stored in a library; from there they can be loaded when needed. The extra code tries to load the package using R’s require() function but\u00a0if that fails, the most likely reason is that the package has never been installed into the user’s R library, so an attempt is made to do that. In my example, I try to install from the R mirror site at Imperial College in London. Good practice says that you ought to install from\u00a0your closest mirror site, so if your intended Stata user is not in southern England then you ought to replace this address with the address of a mirror site in their own country. You can find the addresses of the mirror sites by issuing the following command in R ,<\/p>\n
>\u00a0M <- getCRANmirrors()<\/span><\/p>\nand then looking at M. If your program is intended for Stata users in many different countries then you would need to make the Stata user choose a mirror site, perhaps by creating an ado file that shows the user a list of available sites that you obtain via the getCRANmirrors function.<\/p>\n
I have also written a program called rtemplate.ado<\/strong> that fills in the missing elements of an R template with its arguments. So the command,<\/p>\nrtemplate DPcdensity.txt temp.R a b c \u201c3,4\u201d<\/span><\/p>\nwould read the template DPcdensity.txt from RTEMPLATES and create a file temp.R in the current working directory.\u00a0In temp.R,\u00a0the first @ is replaced by \u201ca\u201d the second by \u201cb\u201d etc. At the same time pounds are converted into dollars. The source for this ado file is very simple,<\/p>\n
program rtemplate<\/span>
\n\u00a0\u00a0\u00a0\u00a0 gettoken template 0 : 0<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>gettoken outfile 0 : 0<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>tempname fhin fhout<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local personal : sysdir PERSONAL<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>file open `fhout’ using “`outfile'”, write replace<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>file open `fhin’ using “`personal’RTEMPLATES\/`template'”, read<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>file read `fhin’ inLine<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>while !r(eof) {<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>tokenize `”`inLine'”‘ , parse(“@\u00a3”)<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>local j = 1<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>while `”“j””‘ != “” {<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>if `”“j””‘ == “@” {<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0 <\/span><\/span>gettoken item 0 : 0<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0 <\/span><\/span>file write `fhout’ `”`item'”‘<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0 <\/span><\/span>}<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0 <\/span><\/span>else if `”“j””‘ == “\u00a3” file write `fhout’ “$”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>else file write `fhout’ `”“j””‘<\/span>
\n \u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0 <\/span><\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>local ++j<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>}<\/span>
\n \u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0 <\/span><\/span>file write `fhout’ _n<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>file read `fhin’ inLine<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>}<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>file close `fhin’<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>file close `fhout’<\/span>
\n end<\/span><\/p>\nNow all I need to do is to write a program that calculates the missing items, @, and sends them to rtemplate.ado in the correct order. For instance,<\/p>\n
program r_DPcdensity<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>syntax varlist [ using\/ ], XPredict(numlist) PMean(name) \/\/\/<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>PCOVariance(name) PAlpha(real) PK0(real) PNU1(real) \/\/\/<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>[ Burnin(integer 1000) Updates(integer 5000) Skip(integer 3) \/\/\/<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>NGrid(integer 50)\u00a0 ] LOGfile(string)<\/span><\/span><\/p>\n* current working directory<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>if “`using'” == “” tempfile using<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local wd = subinstr(“`c(pwd)'”,”\\”,”\/”,.)<\/span>
\n * temporary data file contain x and y<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>tokenize `varlist’<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local y “`1′”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local x “`2′”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>tempfile dfile<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local dfile = subinstr(“`dfile'”,”\\”,”\/”,.)<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>preserve<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>keep `varlist’<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>save “`dfile'”, replace<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>restore<\/span>
\n * mean of base distribution<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local a1 = `pmean'[1,1]<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local a2 = `pmean'[1,2]<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local mv = “`a1′,`a2′”<\/span>
\n * covariance matrix of base distribution<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local a1 = `pcovariance'[1,1]<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local a2 = `pcovariance'[1,2]<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local a3 = `pcovariance'[2,2]<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local cv = “`a1′,`a2′,`a2′,`a3′”<\/span>
\n * df of base distribution<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local df = `pnu1′ – 2 – 1<\/span>
\n * values of xpredict & names for m,l,h<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local xlist “”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local mlist “”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local llist “”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local hlist “”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>local i = 0<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>foreach v of numlist `xpredict’ {<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>local ++i<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>local mlist `”`mlist’,”m`i'””‘<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>local llist `”`llist’,”l`i'””‘<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>local hlist `”`hlist’,”h`i'””‘<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>if “`xlist'” == “” local xlist “`v'”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>else local xlist “`xlist’,`v'”<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>}<\/span>
\n * make the script<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>rtemplate DPcdensity “`using'” “`wd'” “`dfile'” “`mv'” “`cv'” \/\/\/<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>“`palpha'” “`pk0′” “`pnu1′” “`df'” “`burnin'” “`updates'” \/\/\/<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>\u00a0\u00a0\u00a0\u00a0 <\/span>“`skip'” “`xlist'” “`y'” “`x'” “`ngrid'” `”`y’`mlist’`llist’`hlist'”‘<\/span>
\n * run the R script<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>cap rm `logfile’<\/span>
\n \u00a0\u00a0\u00a0\u00a0 <\/span>execute R CMD BATCH `using’ `logfile’<\/span>
\n end<\/span><\/p>\nNow a Stata user can have all of the complexity of R hidden from them and all they need to do is to type commands such as,<\/p>\n
.\u00a0use igg.dta, clear<\/span>
\n. matrix mu = (7,3)<\/span>
\n. matrix sigma = (4,1\\1,1)<\/span>
\n. r_DPcdensity igg age using temp.R , xp(1 2 3 4 5) priorm(mu) priorcov(sigma) \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 priora(5) priork0(0.25) priornu1(10) log(log.out)<\/span>
\n. type log.out<\/span>
\n. use nclusters.dta, clear<\/span>
\n. histogram nclus , start(0.5) width(1) freq title(Number of Clusters)<\/span><\/p>\nand they would fit a Bayesian DPM model with their chosen priors and plot a histogram of the number of clusters at each iteration of the sampler. The computation will be done in R but the Stata user need know nothing of that.<\/p>\n
The key to making the facilities of R available within Stata it to prepare good quality ado files with a structure similar to r_DPcdensity.ado and good quality help files that explain the options. The ado file given above is very basic and would benefit from more checking of the options and some helpful error messages. None the less, writing a good version of r_DPcdensity.ado is very, very much easier than rewriting the source code of the R\u00a0function DPcdensity.<\/p>\n
Next time I will complete this series of postings on using R with Stata by revisiting the way in which results are transferred from R to Stata. Rather than re-package the results into neat .dta files as I have done so far, I will show how the entire set of results can be dumped to a text file and read by Stata.<\/p>\n","protected":false},"excerpt":{"rendered":"
This is another in a series of posting about using 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 […]<\/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":[33,38,4],"class_list":["post-317","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-r","tag-r-packages","tag-stata"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/317","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=317"}],"version-history":[{"count":8,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/317\/revisions"}],"predecessor-version":[{"id":356,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/317\/revisions\/356"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=317"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=317"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=317"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}