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

This is the last 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. I have spent a while on this topic and although I believe it to be very important, next time I want to return to discussing simpler matters.<\/p>\n

Our task for illustrating the process has been to use\u00a0an R package called DPpackage to fit a Bayesian Dirichlet process mixture (DPM) model for smoothing a scatter plot. The dataset for demonstrating the code contains 298 IgG measurements made on children of different ages.<\/p>\n

So far we have seen how to write a Stata ado file that prepares an R script, sends it to R and then reads back the results. To simplify the task, up until now\u00a0I have supposed that the R script includes R commands for packaging the results into Stata .dta files so that reading the results into Stata is simple.<\/p>\n

Many R packages create very complicated list structures of their results and so R code for extraction and re-packaging for export to Stata is not trivial. For this reason I usually adopt a different approach; I dump the entire results structure into a text file and then I use Stata ado files to read the text file and extract any elements that I want. This makes my R scripts simpler and it means that I have access to all of the results in Stata and I am not limited to those exported by the R script.<\/p>\n

First let\u2019s see how we dump the results to a text file. Here is the R template that I developed last time but with the lines that extract the results replaced by new code for dumping the results structure (in green).<\/p>\n

setwd(‘@’)<\/span>
\n\u00a0if( !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><\/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>
\nsink(‘@’)<\/span>
\n dput(fitIgG)<\/span>
\n sink()<\/span><\/p>\n

In R, the sink()\u00a0function redirects output to a file, the name of which will be specified later as it is represented by an @ in the template. Then I write out the\u00a0list structure of results\u00a0using the dput()\u00a0function and finally I redirect the output back to the log file. This approach makes for much simpler R code and places the emphasis of the programming where it belongs, in Stata.<\/p>\n

This\u00a0code saves the entire set of results to a text file, but it is not a text file that is easy to use because the structure of results is quite complex. To illustrate this here are the first few lines of the results created by this particular R script.<\/p>\n

structure(list(call = quote(DPcdensity.default(y = D$igg, x = D$age,<\/span>
\n xpred = xpred, ngrid = 50, compute.band = TRUE, type.band = “PD”,<\/span>
\n prior = prior, mcmc = mcmc, state = NULL, status = TRUE)),<\/span>
\n y = c(1.5, 2.70000004768372, 1.89999997615814, 4, 1.89999997615814,<\/span>
\n 4.40000009536743, 1.5, 2.20000004768372, 1.60000002384186,<\/span>
\n 4.69999980926514, 1.60000002384186, 1.39999997615814, 3,<\/span><\/p>\n

The full text file has 7640 lines.<\/p>\n

To enable me to extract data from such files I have written two Stata ado files. One called rdescribe.ado<\/strong> that reads dumped results and lists all of the items that could be extracted and one called rextract.ado<\/strong> that extracts a requested item.<\/p>\n

Let\u2019s suppose that I dumped my DPcdensity results to a file called fitigg.txt. The Stata command<\/p>\n

.\u00a0rdescribe using fitigg.txt<\/span>
\nwould list the names of the items in the results file but not their values. For this example the list is quite long, so I just give the first few lines.<\/p>\n

Structures and sub-structures within file: fitigg.txt<\/span>
\ncall \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (string)<\/span>
\ny\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single variable)<\/span>
\ncoefficients\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (matrix or set of variables)<\/span>
\nvarnames \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single variable)<\/span>
\nmodelname \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 cpo \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single variable)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 fso \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single variable)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 prior\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(list)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 alpha\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 k0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 m1 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single variable)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 nu1 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 psiinv1 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (matrix or set of variables)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Dim: 2, 2<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 mcmc \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (list)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 nburn\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 nsave\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 nskip\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ndisplay \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 state \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (list)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 alpha \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (single value)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 m1 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (matrix or set of variables)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Dim: 2, 1<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 muclus \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (matrix or set of variables)<\/span>
\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 Dim: 398, 2<\/span><\/p>\n

We see that the results include an item called call<\/strong>, which is a string containing the R command that was executed, next is an item containing the values of the y variable and so it goes on.<\/p>\n

The results structure called prior<\/strong> contains the priors that I chose and is itself a list including items such as alpha<\/strong> and k0<\/strong>. In fact you will see that there are two items called alpha<\/strong>, one in prior<\/strong> and another in state<\/strong> (state gives the contents of the sampler at the end of the run and can be used to continue the chain if desired).<\/p>\n

Once we know what is in the results file we can use rextract.ado<\/strong> to extract any of the items. For instance,
\n\u00a0. rextract using iggresults.txt , item(y) clear<\/span>
\nWill clear any existing data from Stata\u2019s memory and extract y into a Stata variable with that name.<\/p>\n

.\u00a0rextract using iggresults.txt , item(prior(alpha))<\/span>
\nExtracts the value of alpha (a single number) from the list called prior and returns it as the local r(prior_alpha).<\/p>\n

.\u00a0rextract using iggresults.txt , item(state(muclus)) matrix<\/span>
\nExtracts muclus (cluster means at the end of the run) and stores them in a 398×2 Stata matrix called state_mclus.<\/p>\n

.\u00a0rextract using iggresults.txt , item(state(muclus)) clear<\/span>
\nAlso extracts muclus but this time stores the results in two Stata variables called state_muclus_1 and state_muclus_2.<\/p>\n

So there you have it. There is no need to program advanced statistical analyses in Stata or even in Mata. Instead we only\u00a0need to prepare a conversion program that enables a Stata user to select their options and that then prepares the script and sends it to R. In this way any Stata user can have access to these advanced analyses even if they have no idea how to program in R.<\/p>\n

The code for rdescribe.ado<\/strong> and rextract.ado<\/strong> are contained in the pdf rdescribe_rextract<\/a>\u00a0that you can download from here and then copy and paste into separate ado files. The programs were prepared for personal use rather than general release, but you should not have a problem getting them to work.<\/p>\n","protected":false},"excerpt":{"rendered":"

This is the last 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. I have spent a while on this topic and although I believe it to be very important, next time I want to return to […]<\/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,4],"class_list":["post-323","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-r","tag-stata"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/323","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=323"}],"version-history":[{"count":9,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/323\/revisions"}],"predecessor-version":[{"id":363,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/323\/revisions\/363"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=323"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=323"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=323"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}