[R] problem with plots with short example.

[Nicole Ford]

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.