Posted on June 17, 2019 by matsukat

院：認識情報解析　T-test & ANOVA

source("http://www.matsuka.info/univ/course_folder/HDI_revised.txt")
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/TwoGroupIQ.csv")
dat2sd <- dat[dat$Group=="Smart Drug",]

dataList <- list(y = dat2sd$Score, Ntotal = nrow(dat2sd),
                 meanY = mean(dat2sd$Score), sdY = sd(dat2sd$Score))

model.txt = "
model {
  for (i in 1:Ntotal){
    y[i] ~ dnorm(mu, 1/sigma^2)
  }
  mu ~ dnorm(meanY, 1/(100*sdY)^2)
  sigma ~ dunif(sdY/1000, sdY*1000)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma")
library("rjags")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
traceplot(codaSamples)
gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)

HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])
hist(dat2sd$Score, probability = T)
x.temp = seq(40,220,0.1)
n.plot = 100
rand.order = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  y.temp = dnorm(x.temp, mean = mcmcMat[rand.order[i.plot],1],
                 sd = mcmcMat[rand.order[i.plot],2])
  lines(x.temp, y.temp, type='l',col="orange")
}

#y = c(-2:2,15);Ntotal = length(y);meanY = mean(y); sdY = sd(y)
#dataList = list(y= y, Ntotal=Ntotal, meanY = meanY, sdY = sdY)
dataList <- list(y = dat2sd$Score, Ntotal = nrow(dat2sd),
                 meanY = mean(dat2sd$Score), sdY = sd(dat2sd$Score))

model.txt ="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu , 1/sigma^2 , nu )
  }
  mu ~ dnorm( meanY , 1/(100*sdY)^2 )
  sigma ~ dunif( sdY/1000 , sdY*1000 )
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma", "nu")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
#traceplot(codaSamples)
#gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])

dataList <- list(y = dat$Score, Ntotal = nrow(dat),
                 meanY = mean(dat$Score), sdY = sd(dat$Score),
                 Ngroup = length(unique(dat$Group)),
                 group = dat$Group)

model.txt ="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu[group[i]] , 1/sigma[group[i]]^2 , nu )
  }
  for (j in 1:Ngroup){
    mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
    sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
  }
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma", "nu")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
#traceplot(codaSamples)
#gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,1]-mcmcMat[,2])
HDI.plot(mcmcMat[,2])


dat<-read.csv("http://www.matsuka.info/univ/course_folder/FruitflyDataReduced.csv")
y=dat$Longevity
yMean = mean(y)
ySD = sd(y)
Ntotal = length(y)

MeanSD2gamma <- function( mean, sd ) {
  shape = mean^2 / sd^2
  rate = mean / sd^2
  return(data.frame(shape,rate))
}

ModeSD2gamma <- function( mode, sd ) {
  rate = ( mode + sqrt( mode^2 + 4 * sd^2 ) )/( 2 * sd^2 )
  shape = 1 + mode * rate
  return(data.frame(shape,rate))
}

temp.gamma = ModeSD2gamma( mode=sd(y)/2 , sd=2*sd(y) )
gShape = temp.gamma$shape
gRate = temp.gamma$rate

x=as.numeric(dat$CompanionNumber)
Ngroup = length(unique(x))
xc=dat$Thorax
xcMean=mean(xc)

dataList = list(
  y = y ,
  x = x ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  gShape = gShape,
  gRate  = gRate
)

model.txt="
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dnorm( a0 + a[ x[ i_data ] ], 1/ySigma^2)
  }

  ySigma ~ dunif( ySD/100, ySD*10 )
  a0 ~ dnorm( yMean, 1/(ySD*5)^2)

  for ( i_group in 1:Ngroup ) {
    a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
  }
  aSigma ~ dgamma( gShape, gRate )

  for ( i_group in 1:Ngroup ) { m[ i_group ] <- a0 + a[ i_group ] }
  b0 <- mean( m[ 1:Ngroup ] )
  for (i_data in 1:Ngroup ) { b[ i_data ] <- m[ i_data ] - b0 }
  }
"
writeLines(model.txt, "model.txt")

parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma" )

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# fig 19.3
plot(x,dat$Longevity,xlim=c(0,5.1),ylim=c(0,120))
for (i_c in 1:100) {
  sd=mcmcMat[i_c,"ySigma"]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,2:6]+mcmcMat[i_c,"b0"]
  for (i_group in 1:5){
    lower=means[i_group]+lim[1]*sd
    upper=means[i_group]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_group])/sd)
    dens.y=dens.y/max(dens.y)*0.75
    lines(-dens.y+i_group,yseq,type='l',col="orange")
  }
}




dataList = list(
  y = y ,
  x = x ,
  xc = xc ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  xcMean = xcMean ,
  acSD = sd(xc) ,
  gShape = gShape,
  gRate  = gRate
)
model.txt="
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dnorm( mu[i_data], 1/ySigma^2)
    mu [ i_data] <- a0 + a[ x[ i_data ] ] + ac*( xc[ i_data ] - xcMean )
  }

  ySigma ~ dunif( ySD/100, ySD*10 )
  a0 ~ dnorm( yMean, 1/(ySD*5)^2)

  for ( i_group in 1:Ngroup ) {
    a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
  }
  aSigma ~ dgamma( gShape, gRate )
  ac ~ dnorm(0, 1/(2*ySD/acSD)^2 )
  b0 <- a0 + mean( a[ 1:Ngroup ] ) - ac*xcMean
  for (i_data in 1:Ngroup ) { b[ i_data ] <- a[ i_data ] - mean( a[1:Ngroup]) }
}"
writeLines(model.txt, "model.txt")
parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma", "ac" )

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=5000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# fig 19.5 only "none0" will be plotted
plot(Longevity~Thorax,xlim=c(0.6,1),ylim=c(0,120),data=dat[dat$CompanionNumber=="None0",],type='n')
Thorax.value=c(0.65,0.8,0.95)
for (i_c in 1:100) {

  # regression lines
  xs=c(0.5,1.1)
  ys=mcmcMat[i_c,"b0"]+mcmcMat[i_c,3]+mcmcMat[i_c,"ac"]*xs
  lines(xs,ys,col="orange")

  # density
  sd=mcmcMat[i_c,"ySigma"]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,"b0"]+mcmcMat[i_c,3]+mcmcMat[i_c,"ac"]*Thorax.value
  for (i_thorax in 1:3){
    lower=means[i_thorax]+lim[1]*sd
    upper=means[i_thorax]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_thorax])/sd)
    dens.y=dens.y/max(dens.y)*0.05
    lines(-dens.y+Thorax.value[i_thorax],yseq,type='l',col="orange")
  }
}
abline(v=Thorax.value,col='green')
points(Longevity~Thorax,xlim=c(0.6,1),ylim=c(0,120),data=dat[dat$CompanionNumber=="None0",])



dat<-read.csv("http://www.matsuka.info/univ/course_folder/NonhomogVarData.csv")
y=dat$Y
yMean = mean(y)
ySD = sd(y)
Ntotal = length(y)
x=as.numeric(dat$Group)
Ngroup = length(unique(x))

model.txt = "
model {
  for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dt( a0 + a[ x[ i_data ] ], 1/ySigma[x[i_data]]^2, nu)
}
nu <- nuMinusOne + 1
nuMinusOne ~ dexp(1/29)
for (i_group in 1:Ngroup) {
ySigma[i_group]~dgamma(ySigSh,ySigR)
}
ySigSh <- 1+ySigMode * ySigR
ySigR  <- ((ySigMode + sqrt(ySigMode^2+4*ySigSD^2))/(2*ySigSD^2))

ySigMode ~ dgamma(gShape,gRate)
ySigSD ~ dgamma(gShape,gRate)
a0 ~ dnorm( yMean, 1/(ySD*10)^2)

for ( i_group in 1:Ngroup ) {
a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
}
aSigma ~ dgamma( gShape, gRate )

for ( i_group in 1:Ngroup ) {
m[ i_group ] <- a0 + a[ i_group ]
}

b0 <- mean( m[ 1:Ngroup ] )

for (i_data in 1:Ngroup ) {
b[ i_data ] <- m[ i_data ] - b0
}
}"
writeLines(model.txt, "model.txt")

dataList = list(
  y = y ,
  x = x ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  gShape = gShape,
  gRate  = gRate
)

parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma" ,"nu","ySigMode","ySigSD")

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)


# plotting
plot(x,y,xlim=c(0,4.1),ylim=c(70,130))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  i_c = idx[i_sample]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,2:6]+mcmcMat[i_c,"b0"]
  for (i_group in 1:4){
    sd=mcmcMat[i_c,9+i_group]
    lower=means[i_group]+lim[1]*sd
    upper=means[i_group]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_group])/sd)
    dens.y=dens.y/max(dens.y)*0.75
    lines(-dens.y+i_group,yseq,type='l',col="orange")
  }
}

Related