problem with plots with short example.

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130328/d7c69ed2/attachment.pl>
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.

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.
-----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
thank you, dan.  any information, no matter how small, is helpful.

i deleted R this moring and reinstalled it.  outside of that i am not sure what else to delete/ reinstall.  You mention you reinstalled JAGS-  i will give that a try, as well.  thanks!

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

Hi Nicole

I just upgraded to 2.15.3 today
I was just having similar problems with run.jags  and it stopping 
dead in its tracks with an error message pointing to somewhere else

See

?run.jags
and the second paragraph of it.

I then got it to run using a run.jags script using Rterm and saved the model.
By accident I copied and pasted a script to R with run.jags in it and 
it now runs.

see also ?test.jags

Do not know what will happen after a reboot.

HTH

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au
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.

______________________________________________
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.
I was up till 4 am with this sucker trying to figure it out-  I thought i lost my mind/ screwed it up somehow.  I guess it's nice to know it really wasn't me.  Though, I AM sad R is messing up-  first time in over 3 years, so I guess it isn't so bad.  But timing couldn't have been worse, as I have a conference coming up.

Thanks very much, Duncan!  I am going to give that a whirl tomorrow!

Hi Nicole

I just upgraded to 2.15.3 today
I was just having similar problems with run.jags  and it stopping dead in its tracks with an error message pointing to somewhere else

See

?run.jags
and the second paragraph of it.

I then got it to run using a run.jags script using Rterm and saved the model.
By accident I copied and pasted a script to R with run.jags in it and it now runs.

see also ?test.jags

Do not know what will happen after a reboot.

HTH

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au

At 10:24 29/03/2013, you wrote:
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.

______________________________________________
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.
Hi Nicole

My code works using source file if I put

   library(runjags)
   x = testjags(findjags("windows",look_in = "c:/Program Files/JAGS"))
source("G:/Sweave/Bayes.R")

works ok

Further developments with Sweave. Putting this in the Sweave file

   library(runjags)
   x = testjags(findjags("windows",look_in = "c:/Program Files/JAGS"))

x
...
$JAGS.available
[1] TRUE

$JAGS.path
[1] "\"c:/Program Files/JAGS/JAGS-3.3.0/i386/bin/jags-terminal.exe\""
...

If I do not setwd() everything works fine for runjags but have not 
tried to plot any plots as pdf
Unfortunately it comes up in my  default
c:\Users\d mackay\Documents\

############################
Using setwd() before Sweave("G:/Sweave/Bayes.Rnw")
and in the Sweave file

# \Sweave file before run.jags call
   x = testjags(findjags("windows",look_in = "c:/Program Files/JAGS"))

It returns this error in the R command window

Error:  chunk 7 (label = JAGS2)
Error in run.jags(model = mod, monitor = "p", init = list(init1, init2),  :
   Unable to call JAGS

However x shows that it is finding JAGS
x
...
$JAGS.available
[1] TRUE

$JAGS.path
[1] "\"c:/Program Files/JAGS/JAGS-3.3.0/i386/bin/jags-terminal.exe\""

...

Latex tex file
You are currently logged on as d mackay, on a windows machine
You are using R version 2.15.3 (2013-03-01), with the Rgui GUI
JAGS version 3.3.0 found successfully using the command "c:/Program
Files/JAGS/JAGS-3.3.0/i386/bin/jags-terminal.exe"

So do not use setwd or change directory from the gui menu and it 
works but things end up where your default directory is

I am not sure what is going on but any assistance to fix it would be 
appreciated.

Regards

Duncan
I was up till 4 am with this sucker trying to figure it out-  I 
thought i lost my mind/ screwed it up somehow.  I guess it's nice to 
know it really wasn't me.  Though, I AM sad R is messing up-  first 
time in over 3 years, so I guess it isn't so bad.  But timing 
couldn't have been worse, as I have a conference coming up.

Thanks very much, Duncan!  I am going to give that a whirl tomorrow!

On Mar 28, 2013, at 10:37 PM, Duncan Mackay wrote:

Hi Nicole

I just upgraded to 2.15.3 today
I was just having similar problems with run.jags  and it stopping 
dead in its tracks with an error message pointing to somewhere else
See

?run.jags
and the second paragraph of it.

I then got it to run using a run.jags script using Rterm and 
saved the model.
By accident I copied and pasted a script to R with run.jags in it 
and it now runs.
see also ?test.jags

Do not know what will happen after a reboot.

HTH

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au

At 10:24 29/03/2013, you wrote:
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.

______________________________________________
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.
yes I was also having something like this happen to me just before everything blew up. i never set a working directory, but use file.choose() or call the path directly.   but the last few hours before it failed, it kept running older datasets, when i clearly called the new data set, for example,
within.mod <- lm(wvsAB$dem ~ wvsAB$cpi + as.factor(wvsAB$country), data=wvsAB)
would yeild in the summary immediately after:  (with no error called on my model!)
summary(within.mod)
Call:
lm(formula = dat$dem ~ dat$cpi + as.factor(dat$country), 
    data = dat)

i was forced to do rm(data)  more than once.

no idea what's goign on.

Hi Nicole

My code works using source file if I put

 library(runjags)
 x = testjags(findjags("windows",look_in = "c:/Program Files/JAGS"))
source("G:/Sweave/Bayes.R")

works ok

Further developments with Sweave. Putting this in the Sweave file

 library(runjags)
 x = testjags(findjags("windows",look_in = "c:/Program Files/JAGS"))

x
...
$JAGS.available
[1] TRUE

$JAGS.path
[1] "\"c:/Program Files/JAGS/JAGS-3.3.0/i386/bin/jags-terminal.exe\""
...

If I do not setwd() everything works fine for runjags but have not tried to plot any plots as pdf
Unfortunately it comes up in my  default
c:\Users\d mackay\Documents\

############################
Using setwd() before Sweave("G:/Sweave/Bayes.Rnw")
and in the Sweave file

# \Sweave file before run.jags call
 x = testjags(findjags("windows",look_in = "c:/Program Files/JAGS"))

It returns this error in the R command window

Error:  chunk 7 (label = JAGS2)
Error in run.jags(model = mod, monitor = "p", init = list(init1, init2),  :
 Unable to call JAGS

However x shows that it is finding JAGS
x
...
$JAGS.available
[1] TRUE

$JAGS.path
[1] "\"c:/Program Files/JAGS/JAGS-3.3.0/i386/bin/jags-terminal.exe\""

...

Latex tex file
You are currently logged on as d mackay, on a windows machine
You are using R version 2.15.3 (2013-03-01), with the Rgui GUI
JAGS version 3.3.0 found successfully using the command "c:/Program
Files/JAGS/JAGS-3.3.0/i386/bin/jags-terminal.exe"

So do not use setwd or change directory from the gui menu and it works but things end up where your default directory is

I am not sure what is going on but any assistance to fix it would be appreciated.

Regards

Duncan

At 12:45 29/03/2013, you wrote:
I was up till 4 am with this sucker trying to figure it out-  I thought i lost my mind/ screwed it up somehow.  I guess it's nice to know it really wasn't me.  Though, I AM sad R is messing up-  first time in over 3 years, so I guess it isn't so bad.  But timing couldn't have been worse, as I have a conference coming up.

Thanks very much, Duncan!  I am going to give that a whirl tomorrow!

On Mar 28, 2013, at 10:37 PM, Duncan Mackay wrote:

Hi Nicole

I just upgraded to 2.15.3 today
I was just having similar problems with run.jags  and it stopping dead in its tracks with an error message pointing to somewhere else

See

?run.jags
and the second paragraph of it.

I then got it to run using a run.jags script using Rterm and saved the model.
By accident I copied and pasted a script to R with run.jags in it and it now runs.

see also ?test.jags

Do not know what will happen after a reboot.

HTH

Duncan

Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au

At 10:24 29/03/2013, you wrote:
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.