[R] problem with plots with short example.
Nicole Ford
nicole.ford at me.com
Fri Mar 29 03:57:11 CET 2013
LOL dennis, thanks for the laugh at the end.
I'm a PhD student... I had never used it before until this last semester when i started taking Bayes at my new university (transfer because of husband's job =\). My old dept offered bayes, but not the new one, so i am taking it in the psychology dept.... the application (lots of binary outcomes and t-test focused) isn't terribly useful for my own discipline (poli sci, where ALL i do is LM and logistic and now HLM!) but it -at least- has exposed me to how to use bayes, JAGS, etc. so it isn't completely foreign to me.
I am on a mac and, for whatever reason (until the latest update), I was NEVER able to run any code on it. it was always always locking up on me.... everyone else in the class was using it (we use the Kruschke book and he has code, etc and recommnds it) so I was determinded to get it to work.... always ended in HOURS and HOURS of frustration and me just saying *screw it* and going to R and completing the task in 30 minutes. =b
but like i said, since they have updated it (was in Jan of Feb?), i haven't had any issues but, MAN was that a frustrating 5/6 months. it really put a sour taste in my mouth.
but, i will say- thank god for R-studio right now. lol totally saved my sanity. but i am determined to fix R...
On Mar 28, 2013, at 10:45 PM, Dennis Murphy wrote:
> Hi Nicole:
>
> On Thu, Mar 28, 2013 at 7:25 PM, Nicole Ford <nicole.ford at me.com> wrote:
>> thank you, dennis, but that's ok! i ran it all successfully yesterday, as well as my own.... any ideas why it would blow up like this and not plot?
>>
> Sorry, no. Given the tenor of some of your posts today, I wondered if
> your computer had contracted a virus. Progressively weird behavior is
> often a sign, as is a noticeable slowdown in performance.
>
>> what's interesting is: was able to run my own research just fine (JUST NOW finished) in R-studio (YUK!!!!), so something is definitely corrupt in the R environment. I wish I could figure it out! i HATE using R-studio. =\
>
> I have certain issues with it, but I like RStudio for documentation
> writing and I'm about to start package development using RStudio
> projects in conjunction with Git. If you use script files instead of
> typing at the command line or copy/pasting from a text editor, it's
> easy to bounce back and forth between the two; the only thing I miss
> is the 'Run selection' feature in R GUI script shells running under
> Windows. If you're an R package developer, it's definitely the way to
> go, especially if you're using Git and devtools.
>
> Why do you hate RStudio so much? (I'm envisioning a Sanka commercial
> with Robert Young from the 70s where the commercial opens with the
> line: Hey, why so tense? :) [I can guess why, but I'll leave it at
> that.]
>
> To quote the great psychologist Cesar Millan, stay calm and assertive :)
> Dennis
>
>>
>>
>> On Mar 28, 2013, at 10:13 PM, Dennis Murphy wrote:
>>
>>> Hi Nicole:
>>>
>>> I, too, ran the script file without error, but a set of warnings.
>>> Here's the output from the R console:
>>>
>>>> source("Lecture7.R.txt")
>>> Loading required package: coda
>>> Loading required package: lattice
>>> Loading required package: parallel
>>> Calling the simulation... (this may take some time)
>>> Simulation complete. Reading coda files...
>>> Coda files loaded successfully
>>> Abstracting pD ... 50000 valid values
>>> Calculating the Gelman-Rubin statistic for 3 variables....
>>> The Gelman-Rubin statistic is below 1.05 for all parameters
>>> Finished running the simulation
>>> Calling the simulation... (this may take some time)
>>> Simulation complete. Reading coda files...
>>> Coda files loaded successfully
>>> Abstracting pD ... 50000 valid values
>>> Calculating the Gelman-Rubin statistic for 69 variables....
>>> The Gelman-Rubin statistic is below 1.05 for all parameters
>>> Finished running the simulation
>>> Calling the simulation... (this may take some time)
>>> Simulation complete. Reading coda files...
>>> Coda files loaded successfully
>>> Abstracting pD ... 100000 valid values
>>> Calculating the Gelman-Rubin statistic for 73 variables....
>>> The Gelman-Rubin statistic is below 1.05 for all parameters
>>> Finished running the simulation
>>> Package `sm', version 2.2-4.1
>>> Copyright (C) 1997, 2000, 2005, 2007, 2008, A.W.Bowman & A.Azzalini
>>> Type help(sm) for summary information
>>>
>>> Loading required package: Matrix
>>>
>>> Attaching package: 'lme4'
>>>
>>> The following object(s) are masked from 'package:coda':
>>>
>>> HPDinterval
>>>
>>> The following object(s) are masked from 'package:stats':
>>>
>>> AIC, BIC
>>>
>>> Calling the simulation... (this may take some time)
>>> Simulation complete. Reading coda files...
>>> Coda files loaded successfully
>>> Calculating the Gelman-Rubin statistic for 9 variables....
>>> The Gelman-Rubin statistic is below 1.05 for all parameters
>>> Finished running the simulation
>>> Calling the simulation... (this may take some time)
>>> Simulation complete. Reading coda files...
>>> Coda files loaded successfully
>>> Calculating the Gelman-Rubin statistic for 10 variables....
>>> The Gelman-Rubin statistic is below 1.05 for all parameters
>>> Finished running the simulation
>>> There were 14 warnings (use warnings() to see them)
>>>
>>> # Unique warnings:
>>>> warnings()[1]
>>> [1] `Use of the 'monitor.deviance', 'monitor.pd', 'monitor.popt' and
>>> 'monitor.pdi' arguments is deprecated - add the special variables
>>> 'dic', 'deviance', 'pd', 'popt' or 'pd.i' to the 'monitor' argument
>>> instead`
>>> [2] No initial value blocks found and n.chains not specified. 2
>>> chains were used.
>>> NULL
>>>
>>> I saved the R workspace for you, but it's over 1Gb. Do you want it? (I
>>> may have enough space in my Dropbox, but if I don't you may want to
>>> suggest an alternative location.)
>>>
>>> Dennis
>>>
>>>
>>> On Thu, Mar 28, 2013 at 4:55 PM, Nicole Ford <nicole.ford at me.com> wrote:
>>>> i am having problem running my own data. yesterday it was working just fine. today it is not. this is the code i was using as an example to follow. this code ALSO worked just fine yesterday, and is no longer working at all. i suspect it is a problem with either my computer or the software, at this point. if THIS won't even run.... something is wrong.
>>>>
>>>> i can assure you this isn't HW.... i know dave, but i am no longer at UW-M and i have never learned HLMs and i am learning this on my own for my own research.
>>>>
>>>> his code is here, along with data. it is short, quick, etc.
>>>>
>>>> http://www.quantoid.net/936/Lecture7.R
>>>>
>>>> ### R code from vignette source 'Lecture7.Rnw'
>>>>
>>>> ###################################################
>>>> ### code chunk number 1: opts
>>>> ###################################################
>>>> options(useFancyQuotes=F)
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 2: data1
>>>> ###################################################
>>>> library(foreign)
>>>> therms <- na.omit(read.dta("http://quantoid.net/936/2008_difftherm.dta"))
>>>> unstate <- unique(therms[,1])
>>>> therms$numstate <- match(therms$state, unstate)
>>>> library(runjags)
>>>> dat <- dump.format(list(
>>>> N = nrow(therms), J=length(unstate),
>>>> y = therms$difftherm,
>>>> numstate = therms$numstate
>>>> ))
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 3: exchange
>>>> ###################################################
>>>> exchange.mod <- "model{
>>>> for(i in 1:N){
>>>> y[i] ~ dnorm(mu, tau)
>>>> }
>>>> mu ~ dnorm(0,.001)
>>>> tau ~ dgamma(.1,.1)
>>>> }"
>>>> exchange.out <- run.jags(exchange.mod,
>>>> data=dat, burnin=10000, sample=50000,
>>>> thin=5, monitor=c("mu", "tau"),
>>>> monitor.deviance=T, monitor.pd=T,
>>>> silent.jags=T)
>>>>
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 4: exchange
>>>> ###################################################
>>>> FE.mod <- "model{
>>>> for(i in 1:N){
>>>> y[i] ~ dnorm(mu[numstate[i]], tau[numstate[i]])
>>>> }
>>>> for(j in 1:J){
>>>> mu[j] ~ dnorm(0,.001)
>>>> tau[j] ~ dgamma(.1,.1)
>>>> }
>>>> }"
>>>> FE.out <- run.jags(FE.mod,
>>>> data=dat, burnin=10000, sample=50000,
>>>> thin=5, monitor=c("mu", "tau"),
>>>> monitor.deviance=T, monitor.pd=T,
>>>> silent.jags=T)
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 5: exchange
>>>> ###################################################
>>>> hier.mod <- "model{
>>>> for(i in 1:N){
>>>> y[i] ~ dnorm(mu[numstate[i]], tau[numstate[i]])
>>>> }
>>>> for(j in 1:J){
>>>> mu[j] ~ dnorm(theta,nu)
>>>> tau[j] ~ dgamma(a,b)
>>>> }
>>>> theta ~ dnorm(0,.01)
>>>> nu ~ dgamma(.1,.1)
>>>> a ~ dunif(0,1000)
>>>> b ~ dunif(0,1000)
>>>> }"
>>>> hier.out <- run.jags(hier.mod,
>>>> data=dat, burnin=10000, sample=100000,
>>>> thin=10, monitor=c("mu", "tau", "theta", "nu", "a", "b"),
>>>> monitor.deviance=T, monitor.pd=T,
>>>> silent.jags=T)
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 6: sums
>>>> ###################################################
>>>> hier.chains <- combine.mcmc(hier.out$mcmc)
>>>> FE.chains <- combine.mcmc(FE.out$mcmc)
>>>> exchange.chains <- combine.mcmc(exchange.out$mcmc)
>>>>
>>>> mu.bar <- apply(FE.chains[, grep("mu\\[", colnames(FE.chains))], 2, mean)
>>>> mu.bar2 <- apply(hier.chains[, grep("mu\\[", colnames(hier.chains))], 2, mean)
>>>> ns <- aggregate(therms$numstate, list(therms$stateabb), length)
>>>> plot(mu.bar, mu.bar2, cex=sqrt(ns[,2])/3,
>>>> xlab = "FE mu[j]",
>>>> ylab = "Hierarchical mu[j]")
>>>> abline(a=0, b=1)
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 7: dotchart
>>>> ###################################################
>>>> fe.mu <- FE.chains[,grep("mu\\[", colnames(FE.chains))]
>>>> fe.ci <- t(apply(fe.mu, 2, quantile, c(.5,.025,.975)))
>>>> rownames(fe.ci) <- unstate
>>>> fe.ci <- fe.ci[order(fe.ci[,1]), ]
>>>> dotchart(fe.ci[order(fe.ci[,1]),1], lcolor="white", pch=16,
>>>> xlim=range(c(fe.ci)))
>>>> segments(fe.ci[,2], 1:34, fe.ci[,3], 1:34)
>>>> mu.ci <- quantile(exchange.chains[,1], c(.5,.025,.975))
>>>> polygon(x=mu.ci[c(2,3,3,2)],
>>>> y = c(-1,-1,36,36),
>>>> col=rgb(128,128,128,100, maxColorValue=255),
>>>> border=NA)
>>>> abline(v=mu.ci[1], lty=2, lwd=2)
>>>> axis(4, at=1:34, labels=ns[match(rownames(fe.ci), ns[,1]),2],
>>>> cex.axis=.75, las=2)
>>>>
>>>>
>>>> ###################################################
>>>> ### code chunk number 8: femeans
>>>> ###################################################
>>>> library(sm)
>>>> sm.density(mu.bar, model="normal")
>>>>
>>>>
>>>> ############################
>>>>
>>>>
>>>>
>>>>
>>>> [[alternative HTML version deleted]]
>>>>
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>
More information about the R-help
mailing list