model {
for ( i in 1:Ncell ) {
y[i] ~ dpois( lambda[i] )
lambda[i] <- exp( a0 + a1[x1[i]] + a2[x2[i]] + a1a2[x1[i],x2[i]] )
}
a0 ~ dnorm( yLogMean , 1/(yLogSD*2)^2 )
for ( j1 in 1:Nx1Lvl ) { a1[j1] ~ dnorm( 0.0 , 1/a1SD^2 ) }
a1SD ~ dgamma(agammaShRa[1],agammaShRa[2])
for ( j2 in 1:Nx2Lvl ) { a2[j2] ~ dnorm( 0.0 , 1/a2SD^2 ) }
a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
for ( j1 in 1:Nx1Lvl ) {
for ( j2 in 1:Nx2Lvl ) {
a1a2[j1,j2] ~ dnorm( 0.0 , 1/a1a2SD^2 )
}
}
a1a2SD ~ dgamma(agammaShRa[1],agammaShRa[2])
# Convert a0,a1[],a2[],a1a2[,] to sum-to-zero b0,b1[],b2[],b1b2[,] :
for ( j1 in 1:Nx1Lvl ) {
for ( j2 in 1:Nx2Lvl ){
m[j1,j2] <- a0 + a1[j1] + a2[j2] + a1a2[j1,j2] # cell means
}
}
b0 <- mean( m[1:Nx1Lvl,1:Nx2Lvl] )
for ( j1 in 1:Nx1Lvl ) { b1[j1] <- mean( m[j1,1:Nx2Lvl] ) - b0 }
for ( j2 in 1:Nx2Lvl ) { b2[j2] <- mean( m[1:Nx1Lvl,j2] ) - b0 }
for ( j1 in 1:Nx1Lvl ) {
for ( j2 in 1:Nx2Lvl ) {
b1b2[j1,j2] <- m[j1,j2] - ( b0 + b1[j1] + b2[j2] )
}
}
# Compute predicted proportions:
for ( j1 in 1:Nx1Lvl ) {
for ( j2 in 1:Nx2Lvl ) {
expm[j1,j2] <- exp(m[j1,j2])
ppx1x2p[j1,j2] <- expm[j1,j2]/sum(expm[1:Nx1Lvl,1:Nx2Lvl])
}
}
for ( j1 in 1:Nx1Lvl ) { ppx1p[j1] <- sum(ppx1x2p[j1,1:Nx2Lvl]) }
for ( j2 in 1:Nx2Lvl ) { ppx2p[j2] <- sum(ppx1x2p[1:Nx1Lvl,j2]) }
}
Error in parse(text = x, srcfile = src): <text>:1:7: unexpected '{' 1: model { ^
mcmcMat = as.matrix(mcmcCoda)
round(rowSums( mcmcMat[ , grep("ppx1x2p",colnames(mcmcMat)) ] ) ,5)
library(reshape2)
myDF = read.csv( file="../DBDA2Eprograms/HairEyeColor.csv" )
myMat = acast(myDF, Eye ~ Hair, value.var="Count" )
Xsq = chisq.test(myMat)
Xsq$observed
round(Xsq$expected,2)
round(Xsq$residuals,2)
Xsq
Black | Blond | Brown | Red | |
---|---|---|---|---|
Blue | 20 | 94 | 84 | 17 |
Brown | 68 | 7 | 119 | 26 |
Green | 5 | 16 | 29 | 14 |
Hazel | 15 | 10 | 54 | 14 |
Black | Blond | Brown | Red | |
---|---|---|---|---|
Blue | 39.22 | 46.12 | 103.87 | 25.79 |
Brown | 40.14 | 47.20 | 106.28 | 26.39 |
Green | 11.68 | 13.73 | 30.92 | 7.68 |
Hazel | 16.97 | 19.95 | 44.93 | 11.15 |
Black | Blond | Brown | Red | |
---|---|---|---|---|
Blue | -3.07 | 7.05 | -1.95 | -1.73 |
Brown | 4.40 | -5.85 | 1.23 | -0.07 |
Green | -1.95 | 0.61 | -0.35 | 2.28 |
Hazel | -0.48 | -2.23 | 1.35 | 0.85 |
Pearson's Chi-squared test data: myMat X-squared = 138.29, df = 9, p-value < 2.2e-16
myDF = read.csv( file="../DBDA2Eprograms/FourByFourCount.csv" )
countMult = 11 # try 3, 7, and 11
myDF$Count = round( myDF$Count * countMult )
myMat = acast( myDF , A ~ B , value.var="Count")
Xsq = chisq.test(myMat)
Xsq$observed
round(Xsq$expected,2)
round(Xsq$residuals,2)
Xsq
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 22 | 22 | 11 | 11 |
A2 | 22 | 22 | 11 | 11 |
A3 | 11 | 11 | 22 | 22 |
A4 | 11 | 11 | 22 | 22 |
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 16.5 | 16.5 | 16.5 | 16.5 |
A2 | 16.5 | 16.5 | 16.5 | 16.5 |
A3 | 16.5 | 16.5 | 16.5 | 16.5 |
A4 | 16.5 | 16.5 | 16.5 | 16.5 |
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 1.35 | 1.35 | -1.35 | -1.35 |
A2 | 1.35 | 1.35 | -1.35 | -1.35 |
A3 | -1.35 | -1.35 | 1.35 | 1.35 |
A4 | -1.35 | -1.35 | 1.35 | 1.35 |
Pearson's Chi-squared test data: myMat X-squared = 29.333, df = 9, p-value = 0.0005692
myDF = read.csv( file="../DBDA2Eprograms/FourByFourCount.csv" )
countMult = 7 # try 3, 7, and 11
myDF$Count = round( myDF$Count * countMult )
myMat = acast( myDF , A ~ B , value.var="Count")
Xsq = chisq.test(myMat)
Xsq$observed
round(Xsq$expected,2)
round(Xsq$residuals,2)
Xsq
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 14 | 14 | 7 | 7 |
A2 | 14 | 14 | 7 | 7 |
A3 | 7 | 7 | 14 | 14 |
A4 | 7 | 7 | 14 | 14 |
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 10.5 | 10.5 | 10.5 | 10.5 |
A2 | 10.5 | 10.5 | 10.5 | 10.5 |
A3 | 10.5 | 10.5 | 10.5 | 10.5 |
A4 | 10.5 | 10.5 | 10.5 | 10.5 |
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 1.08 | 1.08 | -1.08 | -1.08 |
A2 | 1.08 | 1.08 | -1.08 | -1.08 |
A3 | -1.08 | -1.08 | 1.08 | 1.08 |
A4 | -1.08 | -1.08 | 1.08 | 1.08 |
Pearson's Chi-squared test data: myMat X-squared = 18.667, df = 9, p-value = 0.02818
myDF = read.csv( file="../DBDA2Eprograms/FourByFourCount.csv" )
countMult = 3 # try 3, 7, and 11
myDF$Count = round( myDF$Count * countMult )
myMat = acast( myDF , A ~ B , value.var="Count")
Xsq = chisq.test(myMat)
Xsq$observed
round(Xsq$expected,2)
round(Xsq$residuals,2)
Xsq
Warning message: In chisq.test(myMat): Chi-squared approximation may be incorrect
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 6 | 6 | 3 | 3 |
A2 | 6 | 6 | 3 | 3 |
A3 | 3 | 3 | 6 | 6 |
A4 | 3 | 3 | 6 | 6 |
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 4.5 | 4.5 | 4.5 | 4.5 |
A2 | 4.5 | 4.5 | 4.5 | 4.5 |
A3 | 4.5 | 4.5 | 4.5 | 4.5 |
A4 | 4.5 | 4.5 | 4.5 | 4.5 |
B1 | B2 | B3 | B4 | |
---|---|---|---|---|
A1 | 0.71 | 0.71 | -0.71 | -0.71 |
A2 | 0.71 | 0.71 | -0.71 | -0.71 |
A3 | -0.71 | -0.71 | 0.71 | 0.71 |
A4 | -0.71 | -0.71 | 0.71 | 0.71 |
Pearson's Chi-squared test data: myMat X-squared = 8, df = 9, p-value = 0.5341