thr3ads.net - R help - [R] problem with plots with short example. [Mar 2013]

If this information is useful, please help other people find it:
Share via:

Nicole Ford

2013-Mar-28 23:55 UTC

[R] problem with plots with short example.

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]]

Nicole Ford

2013-Mar-29 00:24 UTC

head link

[R] problem with plots with short example.

To be clear everything "runs" with no error message... the only hint
of a problem is at the end of the code: the plot will not fill out/ it is empty.

if anyone has any idea why something like this might happen, i would greatly
appreciate it...  so i can handle it quickly.

thanks in advance.
On Mar 28, 2013, at 7:55 PM, Nicole Ford 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.

Nordlund, Dan (DSHS/RDA)

2013-Mar-29 00:52 UTC

head link

[R] problem with plots with short example.

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Nicole Ford
> Sent: Thursday, March 28, 2013 4:55 PM
> To: r-help help
> Subject: [R] problem with plots with short example.
> 
> 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")
> 
> 
> ############################
> 
> 
> 
Nicole,

I am not going to be much help, other than to say I just downloaded and
Installed the latest versions of JAGS for Windows, and the rjags and sm
packages.  I am running 64-bit Windows 7 with R-2.15.3.  I cut and pasted the
code from your email, and except for a couple of warnings early on about
deprecated some deprecated parameters, the code seems to have run fine.

So something may have broken your environment, and you may need to do some
reinstallation. Maybe someone else will have some better news for you.

Dan

Daniel J. Nordlund
Washington State Department of Social and Health Services
Planning, Performance, and Accountability
Research and Data Analysis Division
Olympia, WA 98504-5204

Nicole Ford

2013-Mar-29 02:57 UTC

head link

[R] problem with plots with short example.

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.
>>

Possibly Parallel Threads

Search for more apparently analagous threads

R help - Mar 2013 - problem with plots with short example.

[R] problem with plots with short example.

[R] problem with plots with short example.

[R] problem with plots with short example.

[R] problem with plots with short example.

Possibly Parallel Threads