[R] problem with plots with short example.

[Nordlund, Dan (DSHS/RDA)]

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