library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")
dat = read.csv("http://peach.l.chiba-u.ac.jp/course_folder/HtWtData30.csv")
x = dat$height
y = dat$weight
modelS.txt = "
data{
Ntotal <- length(y)
xm <- mean(x)
ym <- mean(y)
xsd <- sd(x)
ysd <- sd(y)
for (i in 1:length(y)){
zx[i] <- (x[i] - xm)/xsd
zy[i] <- (y[i] - ym)/ysd
}
}
model{
for (i in 1:Ntotal){
zy[i] ~ dt(zbeta0 + zbeta1 * zx[i], 1/zsigma^2, nu)
}
zbeta0 ~ dnorm(0, 1/10^2)
zbeta1 ~ dnorm(0, 1/10^2)
zsigma ~ dunif(1.03E-3, 1.03E+3)
nu <- nuMinusOne + 1
nuMinusOne ~ dexp(1/29.0)
#Transfrom to original scale:
beta1 <- zbeta1 * ysd/xsd
beta0 <- zbeta0 * ysd + ym -zbeta1*xm*ysd/xsd
sigma <- zsigma * ysd
}
"
writeLines(modelS.txt, "modelS.txt")
dataList = list(x = x, y = y)
jagsModel =jags.model("modelS.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta1", "beta0", "sigma","zbeta0","zbeta1"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
plot(dat$height,dat$weight,xlab="height",ylab ='weight',pch=20,col="orange",cex=5)
par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1])
plot(mcmcMat[,2],mcmcMat[,1], xlab='B1',ylab="B0",pch=20, col='orange')
plot(mcmcMat[,1],mcmcMat[,2], xlab='B0',ylab="B1",pch=20, col='orange')
HDI.plot(mcmcMat[,2])
par(mfrow=c(2,2))
HDI.plot(mcmcMat[,4])
plot(mcmcMat[,5],mcmcMat[,4], xlab='B1',ylab="B0",pch=20, col='orange')
plot(mcmcMat[,4],mcmcMat[,5], xlab='zB0',ylab="zB1",pch=20, col='orange')
HDI.plot(mcmcMat[,5])
n.mcmc = nrow(mcmcMat)
par(mfrow=c(1,1))
temp.x = c(0,80)
temp.y = mcmcMat[1,1]+mcmcMat[1,2]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight')
idx = sample(1:n.mcmc,100)
for (i.plot in 1:length(idx)){
temp.y = mcmcMat[idx[i.plot],1]+mcmcMat[idx[i.plot],2]*temp.x
lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$height,dat$weight,pch=20,cex=4)
par(mfrow=c(1,1))
n2plot = 100
temp.x = c(0,80)
mean.set = seq(50,80,5)
idx=sample(1:nrow(mcmcMat),n2plot)
plot(x,y,xlim=c(50,80),ylim=c(0,400))
for (i_sample in 1:n2plot) {
temp.y = mcmcMat[idx[i_sample],1]+mcmcMat[idx[i_sample],2]*temp.x
lines(temp.x,temp.y,lwd=2,col="orange")
for (i.means in 1:length(mean.set)){
means = mcmcMat[idx[i_sample],1]+mcmcMat[idx[i_sample],2]*mean.set[i.means]
y.seq = seq(0,400,length.out=1000)
dens.y=dt((y.seq-means)/mcmcMat[idx[i_sample],3],29)
dens.y=dens.y/max(dens.y)
lines(dens.y,y.seq,type='l',col="orange")
lines(-dens.y+mean.set[i.means],y.seq,type='l',col="orange")
}
}
points(x,y,pch=20,cex=3)
model.txt = "
data{
Ntotal <- length(y)
xm <- mean(x)
ym <- mean(y)
xsd <- sd(x)
ysd <- sd(y)
for (i in 1:length(y)){
zx[i] <- (x[i] - xm) / xsd
zy[i] <- (y[i] - ym) / ysd
}
}
model{
for (i in 1:Ntotal){
zy[i] ~ dt( zbeta0[s[i]] + zbeta1[s[i]] * zx[i], 1 / zsigma^2, nu)
}
for (j in 1:Nsubj){
zbeta0[j] ~ dnorm( zbeta0mu, 1/(zbeta0sigma)^2)
zbeta1[j] ~ dnorm( zbeta1mu, 1/(zbeta1sigma)^2)
}
zbeta0mu ~ dnorm(0, 1/(10)^2)
zbeta1mu ~ dnorm(0, 1/(10)^2)
zsigma ~ dunif(1.0E-3, 1.0E+3)
zbeta0sigma ~ dunif(1.0E-3, 1.0E+3)
zbeta1sigma ~ dunif(1.0E-3, 1.0E+3)
nu <- nuMinusOne + 1
nuMinusOne ~ dexp(1/29.0)
for (j in 1:Nsubj){
beta1[j] <- zbeta1[j] * ysd / xsd
beta0[j] <- zbeta0[j] * ysd + ym - zbeta1[j] * xm * ysd / xsd
}
beta1mu <- zbeta1mu * ysd / xsd
beta0mu <- zbeta0mu * ysd + ym -zbeta1mu * xm * ysd / xsd
sigma <- zsigma * ysd
}
"
writeLines(model.txt, "model.txt")
dat = read.csv( file="http://peach.l.chiba-u.ac.jp/course_folder/HierLinRegressData.csv" )
x = dat$X
y = dat$Y
s = dat$Subj
dataList = list(x = x , y = y ,s = s , Nsubj = max(s))
jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta1mu", "beta0mu", "sigma", "beta0","beta1"),
n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,1))
temp.x = c(40,90)
temp.y = mcmcMat[1,"beta0mu"]+mcmcMat[1,"beta1mu"]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight',xlim=c(40,90),ylim=c(40,260))
idx = sample(1:nrow(mcmcMat),100)
for (i.plot in 1:length(idx)){
temp.y = mcmcMat[idx[i.plot],"beta0mu"]+mcmcMat[idx[i.plot],"beta1mu"]*temp.x
lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$X,dat$Y,pch=20,cex=4)
par(mfrow=c(5,5))
temp.x = c(40,90)
for (i.s in 1:25){
temp.y = mcmcMat[1,i.s]+mcmcMat[1,i.s+26]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight',xlim=c(40,90),ylim=c(40,260))
idx = sample(1:nrow(mcmcMat),100)
for (i.plot in 1:length(idx)){
temp.y = mcmcMat[idx[i.plot],i.s]+mcmcMat[idx[i.plot],i.s+26]*temp.x
lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$X[which(dat$Subj==i.s)],dat$Y[which(dat$Subj==i.s)],pch=20,cex=2)
}
# ch 18
dat = read.csv( file="httP://www.matsuka.info/univ/course_folder/Guber1999data.csv" )
par(mfrow=c(1,1))
plot(dat[,c(2,5,8)],pch=20,cex=3)
y = dat$SATT
x = as.matrix(cbind(dat$Spend,dat$PrcntTake))
dataList = list(x = x ,y = y, Nx = dim(x)[2], Ntotal = dim(x)[1])
model.txt = "
data{
ym <- mean(y)
ysd <- sd(y)
for (i in 1:Ntotal){
zy[i] <- (y[i] - ym) / ysd
}
for (j in 1:Nx){
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for (i in 1:Ntotal){
zx[i,j] <- (x[i,j] - xm[j])/xsd[j]
}
}
}
model{
for( i in 1:Ntotal){
zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i, 1:Nx]), 1/zsigma^2, nu)
}
zbeta0 ~ dnorm(0, 1/2^2)
for ( j in 1:Nx){
zbeta[j] ~ dnorm(0, 1/2^2)
}
zsigma ~ dunif(1.0E-5, 1.0E+1)
nu <- nuMinusOne + 1
nuMinusOne ~ dexp(1/29.0)
beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] ) * ysd
beta0 <- zbeta0 * ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx]) * ysd
sigma <- zsigma * ysd
}
"
writeLines(model.txt, "model.txt")
jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta[1]", "beta[2]", "beta0", "sigma"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
plot(as.data.frame(mcmcMat))
par(mfrow=c(1,2))
HDI.plot(mcmcMat[,2])
HDI.plot(mcmcMat[,3])
x = as.matrix(cbind(dat$Spend,dat$PrcntTake,dat$Spend*dat$PrcntTake))
dataList = list(x = x ,y = y, Nx = dim(x)[2], Ntotal = dim(x)[1])
jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta[1]", "beta[2]", "beta[3]", "beta0", "sigma"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,2])
HDI.plot(mcmcMat[,3])
HDI.plot(mcmcMat[,4])
empHDI(mcmcMat[,2])
par(mfrow=c(1,1))
perT = seq(0,80,5)
plot(x=c(perT[1],perT[1]),empHDI(mcmcMat[,2]+mcmcMat[,4]*perT[1]),
xlim=c(-5,85),ylim=c(-15,35),type='l',lwd=3,col="orange",
xlab = c("value on Percent Take"),ylab="Slope on Spend")
for (i.plot in 2:length(perT)){
lines(c(perT[i.plot],perT[i.plot]),empHDI(mcmcMat[,2]+mcmcMat[,4]*perT[i.plot]),
lwd=3,col="orange")
}
abline(h=0,lwd=2,lty=2)
par(mfrow=c(1,1))
perT = seq(0,10,0.5)
plot(x=c(perT[1],perT[1]),empHDI(mcmcMat[,3]+mcmcMat[,4]*perT[1]),
xlim=c(-0.5,11),ylim=c(-5,1),type='l',lwd=3,col="orange",
xlab = c("value on Spent"),ylab="Slope on %Take")
for (i.plot in 2:length(perT)){
lines(c(perT[i.plot],perT[i.plot]),empHDI(mcmcMat[,3]+mcmcMat[,4]*perT[i.plot]),
lwd=3,col="orange")
}
abline(h=0,lwd=2,lty=2)
Related