R/Jags-Ymet-Xnom2grp-MrobustHet.R
In variani/dbda: R code from Doing Bayesian Data Analysis (DBDA)
# Jags-Ymet-Xnom2grp-MrobustHet.R
# Accompanies the book:
# Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
# A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
#source("DBDA2E-utilities.R")
#===============================================================================
genMCMC.JagsYmetXnom2grpMrobustHet = function(object, datFrm, yName="y" , xName="x" ,
saveName=NULL , numSavedSteps=50000 , thinSteps=1 ,
runjagsMethod=runjagsMethodDefault ,
nChains=nChainsDefault ) {
#-----------------------------------------------------------------------------
# THE DATA.
y = as.numeric(datFrm[,yName])
x = as.numeric(as.factor(datFrm[,xName]))
xLevels = levels(as.factor(datFrm[,xName]))
# Do some checking that data make sense:
if ( any( x!=1 & x!=2 ) ) { stop("All x values must be 1 or 2.") }
if ( any( !is.finite(y) ) ) { stop("All y values must be finite.") }
if ( length(x) != length(y) ) { stop("x and y must be same length.") }
Ntotal = length(y)
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
y = y ,
x = x ,
Ntotal = Ntotal ,
meanY = mean(y) ,
sdY = sd(y)
)
#-----------------------------------------------------------------------------
# THE MODEL.
modelString = "
model {
for ( i in 1:Ntotal ) {
y[i] ~ dt( mu[x[i]] , 1/sigma[x[i]]^2 , nu )
}
for ( j in 1:2 ) { # 2 groups
mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
}
nu ~ dexp(1/30.0)
}
" # close quote for modelString
# Write out modelString to a text file
writeLines( modelString , con="TEMPmodel.txt" )
#-----------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
mu = c( mean(y[x==1]) , mean(y[x==2]) )
sigma = c( sd(y[x==1]) , sd(y[x==2]) )
# Regarding initial values in next line: (1) sigma will tend to be too big if
# the data have outliers, and (2) nu starts at 5 as a moderate value. These
# initial values keep the burn-in period moderate.
initsList = list( mu = mu , sigma = sigma , nu = 5 )
#-----------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "mu" , "sigma" , "nu" ) # The parameters to be monitored
adaptSteps = 500 # Number of steps to "tune" the samplers
burnInSteps = 1000
runJagsOut <- run.jags( method=runjagsMethod ,
model="TEMPmodel.txt" ,
monitor=parameters ,
data=dataList ,
inits=initsList ,
n.chains=nChains ,
adapt=adaptSteps ,
burnin=burnInSteps ,
sample=ceiling(numSavedSteps/nChains) ,
thin=thinSteps ,
summarise=FALSE ,
plots=FALSE )
codaSamples = as.mcmc.list( runJagsOut )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
if ( !is.null(saveName) ) {
save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
}
return( codaSamples )
} # end function
#===============================================================================
smryMCMC.JagsYmetXnom2grpMrobustHet = function(object, codaSamples , RopeMuDiff=NULL , RopeSdDiff=NULL ,
RopeEff=NULL , saveName=NULL ) {
summaryInfo = NULL
mcmcMat = as.matrix(codaSamples,chains=TRUE)
summaryInfo = rbind( summaryInfo ,
"mu[1]" = summarizePost( mcmcMat[,"mu[1]"] ) )
summaryInfo = rbind( summaryInfo ,
"mu[2]" = summarizePost( mcmcMat[,"mu[2]"] ) )
summaryInfo = rbind( summaryInfo ,
"muDiff" = summarizePost(
mcmcMat[,"mu[2]"] - mcmcMat[,"mu[1]"] ,
compVal=0.0 , ROPE=RopeMuDiff ) )
summaryInfo = rbind( summaryInfo ,
"sigma[1]" = summarizePost( mcmcMat[,"sigma[1]"] ) )
summaryInfo = rbind( summaryInfo ,
"sigma[2]" = summarizePost( mcmcMat[,"sigma[2]"] ) )
summaryInfo = rbind( summaryInfo ,
"sigmaDiff" = summarizePost(
mcmcMat[,"sigma[2]"] - mcmcMat[,"sigma[1]"] ,
compVal=0.0 , ROPE=RopeSdDiff ) )
summaryInfo = rbind( summaryInfo ,
"nu" = summarizePost( mcmcMat[,"nu"] ) )
summaryInfo = rbind( summaryInfo ,
"log10(nu)" = summarizePost( log10(mcmcMat[,"nu"]) ) )
summaryInfo = rbind( summaryInfo ,
"effSz" = summarizePost(
( mcmcMat[,"mu[2]"]-mcmcMat[,"mu[1]"] )
/ sqrt((mcmcMat[,"sigma[1]"]^2+mcmcMat[,"sigma[2]"]^2)/2) ,
compVal=0.0 , ROPE=RopeEff ) )
if ( !is.null(saveName) ) {
write.csv( summaryInfo , file=paste(saveName,"SummaryInfo.csv",sep="") )
}
return( summaryInfo )
}
#===============================================================================
plotMCMC.JagsYmetXnom2grpMrobustHet = function(object, codaSamples , datFrm , yName="y" , xName="x" ,
RopeMuDiff=NULL , RopeSdDiff=NULL , RopeEff=NULL ,
showCurve=FALSE , pairsPlot=FALSE ,
saveName=NULL , saveType="jpg" ) {
# RopeMuDiff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of means.
# RopeSdDiff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of standard deviations.
# RopeEff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the effect size.
# showCurve is TRUE or FALSE and indicates whether the posterior should
# be displayed as a histogram (by default) or by an approximate curve.
# pairsPlot is TRUE or FALSE and indicates whether scatterplots of pairs
# of parameters should be displayed.
#-----------------------------------------------------------------------------
mcmcMat = as.matrix(codaSamples,chains=TRUE)
chainLength = NROW( mcmcMat )
mu1 = mcmcMat[,"mu[1]"]
mu2 = mcmcMat[,"mu[2]"]
sigma1 = mcmcMat[,"sigma[1]"]
sigma2 = mcmcMat[,"sigma[2]"]
nu = mcmcMat[,"nu"]
#-----------------------------------------------------------------------------
if ( pairsPlot ) {
# Plot the parameters pairwise, to see correlations:
openGraph(width=7,height=7)
nPtToPlot = 1000
plotIdx = floor(seq(1,length(mu1),by=length(mu1)/nPtToPlot))
panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
usr = par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r = (cor(x, y))
txt = format(c(r, 0.123456789), digits=digits)[1]
txt = paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex=1.25 ) # was cex=cex.cor*r
}
pairs( cbind( mu1 , mu2 , sigma1 , sigma2 , log10(nu) )[plotIdx,] ,
labels=c( expression(mu[1]) , expression(mu[2]) ,
expression(sigma[1]) , expression(sigma[2]) ,
expression(log10(nu)) ) ,
lower.panel=panel.cor , col="skyblue" )
if ( !is.null(saveName) ) {
saveGraph( file=paste(saveName,"PostPairs",sep=""), type=saveType)
}
}
#-----------------------------------------------------------------------------
# Set up window and layout:
openGraph(width=6.0,height=8.0)
layout( matrix( c(4,5,7,8,3,1,2,6,9,10) , nrow=5, byrow=FALSE ) )
par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) )
# Select thinned steps in chain for plotting of posterior predictive curves:
nCurvesToPlot = 20
stepIdxVec = seq( 1 , chainLength , floor(chainLength/nCurvesToPlot) )
# Compute limits for plots of data with posterior pred. distributions
y = as.numeric(datFrm[,yName])
x = as.numeric(as.factor(datFrm[,xName]))
xLevels = levels(as.factor(datFrm[,xName]))
y1 = y[x==1]
y2 = y[x==2]
xLim = c( min(y)-0.1*(max(y)-min(y)) , max(y)+0.1*(max(y)-min(y)) )
xBreaks = seq( xLim[1] , xLim[2] ,
length=ceiling((xLim[2]-xLim[1])/(mean(c(sd(y1),sd(y2)))/4)) )
histInfo1 = hist(y1,breaks=xBreaks,plot=FALSE)
histInfo2 = hist(y2,breaks=xBreaks,plot=FALSE)
yMax = 1.2 * max( c( histInfo1$density , histInfo2$density ) )
xVec = seq( xLim[1] , xLim[2] , length=501 )
#-----------------------------------------------------------------------------
# Plot data y1 and smattering of posterior predictive curves:
histInfo = hist( y1 , prob=TRUE , xlim=xLim , ylim=c(0,yMax) , breaks=xBreaks,
col="red2" , border="white" , xlab="y" , ylab="" ,
yaxt="n" , cex.lab=1.5 ,
main=paste("Data for",xLevels[1],"w. Post. Pred.") )
for ( stepIdx in 1:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu1[stepIdxVec[stepIdx]])/sigma1[stepIdxVec[stepIdx]],
df=nu[stepIdxVec[stepIdx]] )/sigma1[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
text( max(xVec) , yMax , bquote(N[1]==.(length(y1))) , adj=c(1.1,1.1) )
#-----------------------------------------------------------------------------
# Plot data y2 and smattering of posterior predictive curves:
histInfo = hist( y2 , prob=TRUE , xlim=xLim , ylim=c(0,yMax) , breaks=xBreaks,
col="red2" , border="white" , xlab="y" , ylab="" ,
yaxt="n" , cex.lab=1.5 ,
main=paste("Data for",xLevels[2],"w. Post. Pred.") )
for ( stepIdx in 1:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu2[stepIdxVec[stepIdx]])/sigma2[stepIdxVec[stepIdx]],
df=nu[stepIdxVec[stepIdx]] )/sigma2[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
text( max(xVec) , yMax , bquote(N[2]==.(length(y2))) , adj=c(1.1,1.1) )
#-----------------------------------------------------------------------------
# Plot posterior distribution of parameter nu:
histInfo = plotPost( log10(nu) , col="skyblue" , # breaks=30 ,
showCurve=showCurve ,
xlab=bquote("log10("*nu*")") , cex.lab = 1.75 ,
cenTend="mode" ,
main="Normality" ) # (<0.7 suggests kurtosis)
#-----------------------------------------------------------------------------
# Plot posterior distribution of parameters mu1, mu2, and their difference:
xlim = range( c( mu1 , mu2 ) )
histInfo = plotPost( mu1 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(mu[1]) , main=paste(xLevels[1],"Mean") ,
col="skyblue" )
histInfo = plotPost( mu2 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(mu[2]) , main=paste(xLevels[2],"Mean") ,
col="skyblue" )
histInfo = plotPost( mu2-mu1 , compVal=0 , showCurve=showCurve ,
xlab=bquote(mu[2] - mu[1]) , cex.lab = 1.75 ,
ROPE=RopeMuDiff,
main="Difference of Means" , col="skyblue" )
#-----------------------------------------------------------------------------
# Plot posterior distribution of param's sigma1, sigma2, and their difference:
xlim=range( c( sigma1 , sigma2 ) )
histInfo = plotPost( sigma1 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(sigma[1]) ,
main=paste(xLevels[1],"Scale") ,
col="skyblue" , cenTend="mode" )
histInfo = plotPost( sigma2 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(sigma[2]) ,
main=paste(xLevels[2],"Scale") ,
col="skyblue" , cenTend="mode" )
histInfo = plotPost( sigma2 - sigma1 ,
compVal=0 , showCurve=showCurve ,
xlab=bquote(sigma[2] - sigma[1]) , cex.lab = 1.75 ,
ROPE=RopeSdDiff ,
main="Difference of Scales" , col="skyblue" ,
cenTend="mode" )
#-----------------------------------------------------------------------------
# Plot effect size.
effectSize = ( mu2 - mu1 ) / sqrt( ( sigma1^2 + sigma2^2 ) / 2 )
histInfo = plotPost( effectSize , compVal=0 , ROPE=RopeEff ,
showCurve=showCurve ,
xlab=bquote( (mu[2]-mu[1])
/sqrt((sigma[1]^2 +sigma[2]^2 )/2 ) ),
cenTend="mode" , cex.lab=1.0 , main="Effect Size" ,
col="skyblue" )
#-----------------------------------------------------------------------------
if ( !is.null(saveName) ) {
saveGraph( file=paste(saveName,"Post",sep=""), type=saveType)
}
}
#===============================================================================
variani/dbda documentation built on May 3, 2019, 4:34 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# Jags-Ymet-Xnom2grp-MrobustHet.R
# Accompanies the book:
# Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:
# A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.
#source("DBDA2E-utilities.R")
#===============================================================================
genMCMC.JagsYmetXnom2grpMrobustHet = function(object, datFrm, yName="y" , xName="x" ,
saveName=NULL , numSavedSteps=50000 , thinSteps=1 ,
runjagsMethod=runjagsMethodDefault ,
nChains=nChainsDefault ) {
#-----------------------------------------------------------------------------
# THE DATA.
y = as.numeric(datFrm[,yName])
x = as.numeric(as.factor(datFrm[,xName]))
xLevels = levels(as.factor(datFrm[,xName]))
# Do some checking that data make sense:
if ( any( x!=1 & x!=2 ) ) { stop("All x values must be 1 or 2.") }
if ( any( !is.finite(y) ) ) { stop("All y values must be finite.") }
if ( length(x) != length(y) ) { stop("x and y must be same length.") }
Ntotal = length(y)
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
y = y ,
x = x ,
Ntotal = Ntotal ,
meanY = mean(y) ,
sdY = sd(y)
)
#-----------------------------------------------------------------------------
# THE MODEL.
modelString = "
model {
for ( i in 1:Ntotal ) {
y[i] ~ dt( mu[x[i]] , 1/sigma[x[i]]^2 , nu )
}
for ( j in 1:2 ) { # 2 groups
mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
}
nu ~ dexp(1/30.0)
}
" # close quote for modelString
# Write out modelString to a text file
writeLines( modelString , con="TEMPmodel.txt" )
#-----------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
mu = c( mean(y[x==1]) , mean(y[x==2]) )
sigma = c( sd(y[x==1]) , sd(y[x==2]) )
# Regarding initial values in next line: (1) sigma will tend to be too big if
# the data have outliers, and (2) nu starts at 5 as a moderate value. These
# initial values keep the burn-in period moderate.
initsList = list( mu = mu , sigma = sigma , nu = 5 )
#-----------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "mu" , "sigma" , "nu" ) # The parameters to be monitored
adaptSteps = 500 # Number of steps to "tune" the samplers
burnInSteps = 1000
runJagsOut <- run.jags( method=runjagsMethod ,
model="TEMPmodel.txt" ,
monitor=parameters ,
data=dataList ,
inits=initsList ,
n.chains=nChains ,
adapt=adaptSteps ,
burnin=burnInSteps ,
sample=ceiling(numSavedSteps/nChains) ,
thin=thinSteps ,
summarise=FALSE ,
plots=FALSE )
codaSamples = as.mcmc.list( runJagsOut )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
if ( !is.null(saveName) ) {
save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
}
return( codaSamples )
} # end function
#===============================================================================
smryMCMC.JagsYmetXnom2grpMrobustHet = function(object, codaSamples , RopeMuDiff=NULL , RopeSdDiff=NULL ,
RopeEff=NULL , saveName=NULL ) {
summaryInfo = NULL
mcmcMat = as.matrix(codaSamples,chains=TRUE)
summaryInfo = rbind( summaryInfo ,
"mu[1]" = summarizePost( mcmcMat[,"mu[1]"] ) )
summaryInfo = rbind( summaryInfo ,
"mu[2]" = summarizePost( mcmcMat[,"mu[2]"] ) )
summaryInfo = rbind( summaryInfo ,
"muDiff" = summarizePost(
mcmcMat[,"mu[2]"] - mcmcMat[,"mu[1]"] ,
compVal=0.0 , ROPE=RopeMuDiff ) )
summaryInfo = rbind( summaryInfo ,
"sigma[1]" = summarizePost( mcmcMat[,"sigma[1]"] ) )
summaryInfo = rbind( summaryInfo ,
"sigma[2]" = summarizePost( mcmcMat[,"sigma[2]"] ) )
summaryInfo = rbind( summaryInfo ,
"sigmaDiff" = summarizePost(
mcmcMat[,"sigma[2]"] - mcmcMat[,"sigma[1]"] ,
compVal=0.0 , ROPE=RopeSdDiff ) )
summaryInfo = rbind( summaryInfo ,
"nu" = summarizePost( mcmcMat[,"nu"] ) )
summaryInfo = rbind( summaryInfo ,
"log10(nu)" = summarizePost( log10(mcmcMat[,"nu"]) ) )
summaryInfo = rbind( summaryInfo ,
"effSz" = summarizePost(
( mcmcMat[,"mu[2]"]-mcmcMat[,"mu[1]"] )
/ sqrt((mcmcMat[,"sigma[1]"]^2+mcmcMat[,"sigma[2]"]^2)/2) ,
compVal=0.0 , ROPE=RopeEff ) )
if ( !is.null(saveName) ) {
write.csv( summaryInfo , file=paste(saveName,"SummaryInfo.csv",sep="") )
}
return( summaryInfo )
}
#===============================================================================
plotMCMC.JagsYmetXnom2grpMrobustHet = function(object, codaSamples , datFrm , yName="y" , xName="x" ,
RopeMuDiff=NULL , RopeSdDiff=NULL , RopeEff=NULL ,
showCurve=FALSE , pairsPlot=FALSE ,
saveName=NULL , saveType="jpg" ) {
# RopeMuDiff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of means.
# RopeSdDiff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the difference of standard deviations.
# RopeEff is a two element vector, such as c(-1,1), specifying the limit
# of the ROPE on the effect size.
# showCurve is TRUE or FALSE and indicates whether the posterior should
# be displayed as a histogram (by default) or by an approximate curve.
# pairsPlot is TRUE or FALSE and indicates whether scatterplots of pairs
# of parameters should be displayed.
#-----------------------------------------------------------------------------
mcmcMat = as.matrix(codaSamples,chains=TRUE)
chainLength = NROW( mcmcMat )
mu1 = mcmcMat[,"mu[1]"]
mu2 = mcmcMat[,"mu[2]"]
sigma1 = mcmcMat[,"sigma[1]"]
sigma2 = mcmcMat[,"sigma[2]"]
nu = mcmcMat[,"nu"]
#-----------------------------------------------------------------------------
if ( pairsPlot ) {
# Plot the parameters pairwise, to see correlations:
openGraph(width=7,height=7)
nPtToPlot = 1000
plotIdx = floor(seq(1,length(mu1),by=length(mu1)/nPtToPlot))
panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
usr = par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r = (cor(x, y))
txt = format(c(r, 0.123456789), digits=digits)[1]
txt = paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex=1.25 ) # was cex=cex.cor*r
}
pairs( cbind( mu1 , mu2 , sigma1 , sigma2 , log10(nu) )[plotIdx,] ,
labels=c( expression(mu[1]) , expression(mu[2]) ,
expression(sigma[1]) , expression(sigma[2]) ,
expression(log10(nu)) ) ,
lower.panel=panel.cor , col="skyblue" )
if ( !is.null(saveName) ) {
saveGraph( file=paste(saveName,"PostPairs",sep=""), type=saveType)
}
}
#-----------------------------------------------------------------------------
# Set up window and layout:
openGraph(width=6.0,height=8.0)
layout( matrix( c(4,5,7,8,3,1,2,6,9,10) , nrow=5, byrow=FALSE ) )
par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) )
# Select thinned steps in chain for plotting of posterior predictive curves:
nCurvesToPlot = 20
stepIdxVec = seq( 1 , chainLength , floor(chainLength/nCurvesToPlot) )
# Compute limits for plots of data with posterior pred. distributions
y = as.numeric(datFrm[,yName])
x = as.numeric(as.factor(datFrm[,xName]))
xLevels = levels(as.factor(datFrm[,xName]))
y1 = y[x==1]
y2 = y[x==2]
xLim = c( min(y)-0.1*(max(y)-min(y)) , max(y)+0.1*(max(y)-min(y)) )
xBreaks = seq( xLim[1] , xLim[2] ,
length=ceiling((xLim[2]-xLim[1])/(mean(c(sd(y1),sd(y2)))/4)) )
histInfo1 = hist(y1,breaks=xBreaks,plot=FALSE)
histInfo2 = hist(y2,breaks=xBreaks,plot=FALSE)
yMax = 1.2 * max( c( histInfo1$density , histInfo2$density ) )
xVec = seq( xLim[1] , xLim[2] , length=501 )
#-----------------------------------------------------------------------------
# Plot data y1 and smattering of posterior predictive curves:
histInfo = hist( y1 , prob=TRUE , xlim=xLim , ylim=c(0,yMax) , breaks=xBreaks,
col="red2" , border="white" , xlab="y" , ylab="" ,
yaxt="n" , cex.lab=1.5 ,
main=paste("Data for",xLevels[1],"w. Post. Pred.") )
for ( stepIdx in 1:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu1[stepIdxVec[stepIdx]])/sigma1[stepIdxVec[stepIdx]],
df=nu[stepIdxVec[stepIdx]] )/sigma1[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
text( max(xVec) , yMax , bquote(N[1]==.(length(y1))) , adj=c(1.1,1.1) )
#-----------------------------------------------------------------------------
# Plot data y2 and smattering of posterior predictive curves:
histInfo = hist( y2 , prob=TRUE , xlim=xLim , ylim=c(0,yMax) , breaks=xBreaks,
col="red2" , border="white" , xlab="y" , ylab="" ,
yaxt="n" , cex.lab=1.5 ,
main=paste("Data for",xLevels[2],"w. Post. Pred.") )
for ( stepIdx in 1:length(stepIdxVec) ) {
lines(xVec, dt( (xVec-mu2[stepIdxVec[stepIdx]])/sigma2[stepIdxVec[stepIdx]],
df=nu[stepIdxVec[stepIdx]] )/sigma2[stepIdxVec[stepIdx]] ,
type="l" , col="skyblue" , lwd=1 )
}
text( max(xVec) , yMax , bquote(N[2]==.(length(y2))) , adj=c(1.1,1.1) )
#-----------------------------------------------------------------------------
# Plot posterior distribution of parameter nu:
histInfo = plotPost( log10(nu) , col="skyblue" , # breaks=30 ,
showCurve=showCurve ,
xlab=bquote("log10("*nu*")") , cex.lab = 1.75 ,
cenTend="mode" ,
main="Normality" ) # (<0.7 suggests kurtosis)
#-----------------------------------------------------------------------------
# Plot posterior distribution of parameters mu1, mu2, and their difference:
xlim = range( c( mu1 , mu2 ) )
histInfo = plotPost( mu1 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(mu[1]) , main=paste(xLevels[1],"Mean") ,
col="skyblue" )
histInfo = plotPost( mu2 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(mu[2]) , main=paste(xLevels[2],"Mean") ,
col="skyblue" )
histInfo = plotPost( mu2-mu1 , compVal=0 , showCurve=showCurve ,
xlab=bquote(mu[2] - mu[1]) , cex.lab = 1.75 ,
ROPE=RopeMuDiff,
main="Difference of Means" , col="skyblue" )
#-----------------------------------------------------------------------------
# Plot posterior distribution of param's sigma1, sigma2, and their difference:
xlim=range( c( sigma1 , sigma2 ) )
histInfo = plotPost( sigma1 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(sigma[1]) ,
main=paste(xLevels[1],"Scale") ,
col="skyblue" , cenTend="mode" )
histInfo = plotPost( sigma2 , xlim=xlim , cex.lab = 1.75 ,
showCurve=showCurve ,
xlab=bquote(sigma[2]) ,
main=paste(xLevels[2],"Scale") ,
col="skyblue" , cenTend="mode" )
histInfo = plotPost( sigma2 - sigma1 ,
compVal=0 , showCurve=showCurve ,
xlab=bquote(sigma[2] - sigma[1]) , cex.lab = 1.75 ,
ROPE=RopeSdDiff ,
main="Difference of Scales" , col="skyblue" ,
cenTend="mode" )
#-----------------------------------------------------------------------------
# Plot effect size.
effectSize = ( mu2 - mu1 ) / sqrt( ( sigma1^2 + sigma2^2 ) / 2 )
histInfo = plotPost( effectSize , compVal=0 , ROPE=RopeEff ,
showCurve=showCurve ,
xlab=bquote( (mu[2]-mu[1])
/sqrt((sigma[1]^2 +sigma[2]^2 )/2 ) ),
cenTend="mode" , cex.lab=1.0 , main="Effect Size" ,
col="skyblue" )
#-----------------------------------------------------------------------------
if ( !is.null(saveName) ) {
saveGraph( file=paste(saveName,"Post",sep=""), type=saveType)
}
}
#===============================================================================
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.