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