There are two ways of implementing a Stata-R link; if you know the R language then you can code the R syntax directly into a commented-out section within your do file and export it to R in just the same way that we export a WinBUGS model file; on the other hand, if you want the R package to be available to users who only know Stata, then you have to write a Stata ado file that writes the R code for them.<\/p>\n
It is instructive to start with a program that assumes some knowledge of R;\u00a0next week\u00a0we will see how we can hide the complexity. I will assume that you have R installed on your computer and that you have installed the package MCMCglmm by downloading it over the internet. This installation\u00a0will also be hidden from the user in next week’s version.<\/p>\n
I’ve discussed using a Stata-R link before in this blog and I am going to follow the same approach. You might find it helpful to read my original postings on this topic as they cover some of the basics in more detail. You can find them at, The data for our Poisson regression model describe the numbers of seizures experienced by 59 epilepsy patients in each of 4 consecutive 2 weeks periods. There are\u00a0five explanatory variables and two random effects; one with a different level for each subject and one with a different level for each subject in each week. The\u00a0problem was\u00a0first introduced to this blog back in February https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2015\/02\/06\/poisson-regression-with-two-random-effects\/<\/a><\/p>\n MCMCglmm<\/strong> wants the data stacked into long columns with the four visits stored under one another. We could re-organize the data in R (R has a very neat package called dplyr<\/strong> that is designed specifically for data organization) but as we are primarily using Stata, we will prepare the data in Stata so that it is already in the format that MCMCglmm needs.<\/p>\n Here is the Stata code that centres the data, stacks it and saves it. In fact, MCMCglmm does not require us to centre the variables but, if I did not, then I would have to change the priors that I used in previous postings, so I have decided to centre for consistency.<\/p>\n use epilepsy.dta , clear<\/span> gen id = _n<\/span> Notice that the data are stored in Stata 12\u00a0format, R cannot read Stata 13 or 14. Now all we need to do is to read these data\u00a0into R and run MCMCglmm. Here is\u00a0the R code that we need. To help avoid confusion, I will show R commands and output in green and Stata code and its output in blue.<\/p>\n library(foreign)<\/span> The library()<\/strong> command (require() is an alternative that does much the same thing) loads two packages; foreign<\/strong>, which reads data files created by other programs and MCMCglmm<\/strong>, which provides our analysis.<\/p>\n This R code\u00a0reads the Stata data file and stores the data under the name D<\/strong>.\u00a0 Then\u00a0it defines the prior and runs the analysis storing the results under the name fit<\/strong>. The summary command prints a table of the estimates. In this case we get<\/p>\n The estimates are similar to those given by WinBUGS.<\/p>\n Our next\u00a0task is to understand the call to MCMCglmm, which is straightforward, and to understand the way that the priors are defined, which is not so straightforward.<\/p>\n The call to MCMCglmm almost explains itself. The fixed part of the model comes first and this is followed by the random effects; in this case I have a single random effect for each\u00a0subject and the structure, R, defined in the prior provides the second random effect with a different value for each observation. The Poisson family\u00a0is given\u00a0a log-link by default. The analysis consists of 55000 iterations but the first 5000 of these are discarded and the remainder are thinned by only saving every 5th simulation.<\/p>\n The first thing to know when trying to understand the priors is that in MCMCglmm all families are\u00a0given a residual variance structure, R, which in this case acts in addition to the randomness of the Poisson distribution. This residual variance is analogous to defining a random effect that takes a different value for each observation.<\/p>\n The other random effects are collectively called G and if there are many of them they are distinguished as G1, G2 \u2026 In this case we only define a random effect for each subject.<\/p>\n The priors on the variances of the random effects are inverse-Wishart distributions; one of the limitations of MCMCglmm\u00a0is that it insists on this distribution as it is conjugate and simplifies the MCMC algorithm.\u00a0When the random effects are independent, as they are here, they\u00a0have a single variance and the multi-dimensional inverse-Wishart distribution reduces to an inverse-gamma. We need to understand this parameterization in order to set the priors.<\/p>\n MCMCglmm requires the user to parameterize the inverse-gamma prior by providing a guess at the true variance, V<\/strong>, and a parameter, nu<\/strong>, that expresses our certainty in the guess, a large value for nu and we are expressing strong belief in the accuracy of our guess at V.<\/p>\n The standard deviation of log(variance) is approximately sqrt(2\/( nu-4)) so if we guess that V=3 and set nu=5 then we have 95% belief that the log(variance) will lie in the range log(3)\u00b12sqrt(2) or equivalently that the variance will lie in the range exp(log(3)\u00b12sqrt(2)) or between 0.2 and 50. Alternatively we could set nu=50 to reflect belief that the variance is most likely to be between 2 and 4.5.<\/p>\n Previously we worked with the Gamma distribution for the precision and set \u03b1=8 and \u03b2=1\/2 so the expected precision was \u03b1\u03b2=4 corresponding to a variance of \u00bc and a standard deviation of \u00bd. The equivalent inverse- gamma just inverts the \u03b2 parameter giving \u03b1=8 and \u03b2=2. The IG has expected variance is \u03b2\/(\u03b1-1)=2\/7 corresponding to a standard deviation of 0.52 (remember that the average(variance) is only approximately 1\/average(precision)).<\/p>\n In the example above I set V=2\/7 and nu=16 to reflect belief that the variance lies between 0.13 and 0.65;\u00a0this is quite informative but similar to my previous analyses of these data. This prior says that I believe that the standard deviation for the random effect will be in the range (0.36,0.80)<\/p>\n Next, we ought to make a tour of the presentation of these results, which I will not try to defend.<\/p>\n The two random variances (0.2607,0.1751) correspond to standard\u00a0deviations (0.51, 0.42); very similar to those that we obtained by WinBUGS.<\/p>\n The idea of quoting 5 decimal places in the table of coefficients is fanciful, probably 2 decimal places are just about reliable given the number of iterations.<\/p>\n The effective sample size is another way of presenting the standard error of the mean. Here the mixing is such that a sample of 10,000 values after thinning acts like a random sample of 5071 when we estimate the coefficient of v4. The standard deviation of the simulated v4 coefficients was 0.09285 so had mixing been perfect so that the simulations\u00a0were independent, the\u00a0standard error would have been 0.09285\/sqrt(10000)=0.0009285 but actually when we allow for the correlation it is 0.001285=0.09285\/sqrt(n) where n is 5071.<\/p>\n The disadvantage of this presentation is that the sem is not easily derived so we cannot tell whether to trust 1,2,3,4 or 5 decimal places in the posterior means.<\/p>\n pMCMC is a complete cheat. It is not a p-value at all despite the stars that R uses to indicate significance, but rather it is twice the probability that the true value is above or below zero. So, in\u00a0the 10,000 simulations v4 must have been negative most of the time but positive on 1465 iterations,\u00a0giving a 0.1465 chance that it is +ve and pMCMC=2*0.1465=0.293. pMCMC does express our confidence in the sign of the coefficient but the presentation invites confusion.<\/p>\n The CI’s are not confidence intervals or even credible intervals but highest posterior density (HPD) intervals.<\/p>\n I’ve explained the output so that we can see the parallels between this analysis and the analyses that we obtained from WinBUGS\u00a0and similar programs\u00a0in previous postings. However, a Stata user need not rely on this presentation, they will read the raw results into Stata and they can summarise them in any way they choose.<\/p>\n The raw results are stored in the R list structure called\u00a0fit<\/strong>. We could write this structure to a text file and then read all or part of that file into Stata. In my previous postings I presented two Stata ado files called rdescribe.ado<\/strong> and rextract.ado<\/strong> that enable us to read R lists from files. The full process comes in two stages first we must add a couple of lines to our R code in order to write fit<\/strong> to a file, then we use rdescribe and rextract to read it in Stata.<\/p>\n The extra R code is,<\/p>\n sink(\u201cmyout.txt\u201d)<\/span> The first call to the function sink() redirects output away from the screen and into a file and the second call returns output to the screen. dput() writes the contents of the structure fit. In fact, we will run R in the background so the default output goes to a log file and not to the screen.<\/p>\n To write the R script to a file and run it we can add the following lines to our Stata do file.<\/p>\n wbsmodel MCMCglmm.do MCMCglmm.R<\/span> This code\u00a0uses\u00a0the wbsmodel<\/strong> command that I wrote for creating WinBUGS model files to copy the R commands to a text file called MCMCglmm.R, then the commands\u00a0are sent to R using\u00a0my execute<\/strong> command. Instructions for downloading all of my commands are given at the end of this posting.<\/p>\n execute<\/strong> assumes that you have set up a file called executables.txt<\/strong> in your personal directory containing a line such as, The raw numerical results are in myout.txt<\/strong>. R places a lot of information in the file and it is formatted as a complex list structure. To see what we have, give the Stata command, Structures and sub-structures within file: myout.txt<\/span> The key information is held in Sol<\/strong> (presumably solution). This contains the 10000 simulations of the 6 regression coefficients. The variances of the random effects are in VCV<\/strong>.<\/p>\n We can read Sol into Stata, rename the variables and plot them using the Stata commands<\/p>\n rextract using myout.txt, item(Sol) clear<\/span><\/p>\n rename Sol_1 cons<\/span> mcmctrace cons base trt bxt age v4<\/span><\/p>\n The mixing is excellent and in log.txt we will find that the timing was good also at about 18 seconds. Better than WinBUGS etc. although of course this program only fits glmms.<\/p>\n
\nhttps:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/09\/19\/using-r-with-stata-part-i\/<\/a>
\nhttps:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/09\/26\/using-r-with-stata-part-ii\/<\/a>
\nhttps:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/10\/03\/using-r-with-stata-part-iii\/<\/a>
\nhttps:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/10\/10\/using-r-with-stata-part-iv\/<\/a><\/p>\n
\n gen LB = log(Base\/4)<\/span>
\n qui su LB<\/span>
\n gen cLB = LB – r(mean)<\/span>
\n gen lnAge = log(Age)<\/span>
\n qui su lnAge<\/span>
\n gen clnAge = lnAge – r(mean)<\/span>
\n gen BXT = LB*Trt<\/span>
\n qui su BXT<\/span>
\n qui gen cBXT = BXT – r(mean)<\/span>
\n qui su Trt<\/span>
\n gen cTrt = Trt – r(mean)<\/span><\/p>\n
\n gen zero = 0<\/span>
\n gen one = 1<\/span>
\n stack id y_1 cLB clnAge cTrt cBXT zero id y_2 cLB clnAge cTrt cBXT zero \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0 id y_3 cLB clnAge cTrt cBXT zero id y_4 cLB clnAge cTrt cBXT one , \/\/\/<\/span>
\n\u00a0\u00a0\u00a0\u00a0 into(subject y base age trt bxt v4)<\/span>
\n rename _stack visit<\/span>
\n sort visit subject<\/span>
\n save epilepsyforR.dta , replace version(12)<\/span><\/p>\n
\n library(MCMCglmm)<\/span>
\n D <- read.dta(file=”K\/epilepsy \/epilepsyforR.dta”)<\/span>
\n prior <- list( B=list(mu=c(2.6,1,0,0,0,0),V=diag(1,6,6)),
\nG=list( G1=list(V=2\/7,n=16) ),R=list(V=2\/7,n=16) )<\/span>
\n fit <- MCMCglmm(y ~ base+trt+bxt+age+v4,random= ~subject,<\/span>
\n \u00a0\u00a0\u00a0\u00a0 prior=prior,family=”poisson”,nitt=55000,burnin=5000,thin=5,data=D)<\/span>
\n summary(fit)<\/span><\/p>\nDIC: 1157.116<\/span><\/strong><\/pre>\nG-structure: ~subject<\/span><\/strong><\/pre>\npost.mean l-95% CI u-95% CI eff.samp\r\n subject 0.2607 0.1566 0.3778 4693<\/span><\/strong><\/pre>\nR-structure: ~units<\/span><\/strong><\/pre>\npost.mean l-95% CI u-95% CI eff.samp\r\n units 0.1751 0.1189 0.2356 3639<\/span><\/strong><\/pre>\nLocation effects: y ~ base + trt + bxt + age + v4<\/span><\/strong><\/pre>\n post.mean l-95% CI u-95% CI eff.samp pMCMC\r\n (Intercept) 1.59213 1.42740 1.75384 5924 <1e-04 ***\r\n base 0.91739 0.64866 1.17659 8145 <1e-04 ***\r\n trt -0.76965 -1.52924 -0.01956 5279 0.046 *\r\n bxt 0.25788 -0.12209 0.64963 6500 0.196\r\n age 0.39816 -0.27944 1.07435 4970 0.256\r\n v4 -0.09743 -0.27435 0.08035 5071 0.293<\/span><\/strong><\/pre>\n
\n dput(fit)<\/span>
\n sink()<\/span><\/p>\n
\n \/*<\/span>
\n library(foreign)<\/span>
\n library(MCMCglmm)<\/span>
\n D <- read.dta(file=”E:\/blog\/topics\/epilepsy\/epilepsyforR.dta”)<\/span>
\n prior <- list( B=list(mu=c(2.6,1,0,0,0,0),V=diag(1,6,6)),<\/span>
\n G=list( G1=list(V=2\/7,n=16) ),R=list(V=2\/7,n=16) )<\/span>
\n fit <- MCMCglmm(y ~ base+trt+bxt+age+v4,random= ~subject,<\/span>
\n prior=prior,family=”poisson”,nitt=55000,burnin=5000,thin=5,data=D)<\/span>
\n summary(fit)<\/span>
\n sink(“myout.txt”)<\/span>
\n dput(fit)<\/span>
\n sink()<\/span>
\n *\/<\/span>
\n execute R CMD BATCH MCMCglmm.R log.txt<\/span><\/p>\n
\nR,”Z:\\Software\\R-3.1.2\\bin\\x64\\R.exe”<\/span>
\nThat is, it points to the location of the R.exe file on your system. The execute<\/strong> command sends the R output to log.txt. So in the log file you will find any error messages or, if everything works, you will find the results as produced by R’s summary(fit) command. The most convenient way to inspect the file is with the Stata command,
\nview log.txt<\/span><\/p>\n
\nrdescribe using myout.txt<\/span>
\nand you will be told that we have,<\/p>\n
\n Sol (matrix or set of variables)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 Dim: 10000, 6<\/span>
\nLambda (single value)<\/span>
\n VCV (matrix or set of variables)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 Dim: 10000, 2<\/span>
\n CP (single value)<\/span>
\n Liab (single value)<\/span>
\n Fixed (list)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 formula (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nfl (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nll (single value)<\/span>
\n Random (list)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 formula (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nfl (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nrl (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nrt (single value)<\/span>
\n Residual (list)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 formula (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nfl (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nrl (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 nrt (single value)<\/span>
\n Deviance (matrix or set of variables)<\/span>
\n DIC (single value)<\/span>
\nX (single value)<\/span>
\n Z (single value)<\/span>
\n XL (single value)<\/span>
\n error.term (single variable)<\/span>
\n family (single variable)<\/span>
\n Tune (single value)<\/span>
\n\u00a0\u00a0\u00a0\u00a0 Dim: 1, 1<\/span><\/p>\n
\n rename Sol_2 base<\/span>
\n rename Sol_3 trt<\/span>
\n rename Sol_4 bxt<\/span>
\n rename Sol_5 age<\/span>
\n rename Sol_6 v4<\/span><\/p>\n
![\"Rtrace\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2015\/04\/Rtrace-e1429007987705.png\")
<\/a><\/p>\n