S-PLUS CODES
USED IN THE BOOK:
STATISTICAL METHODS IN FOOD AND CONSUMER RESEARCH
By M. Gacula, J. Singh, J. Bi, S. Altan
Sensometrics Research and Service
9212 Groomfield Road
Richmond, Virginia 23236, USA
All rights reserved
LIST OF S-PLUS CODES (CHAPTERS 9-11)
No. / Name / Section / Page in the book / Page in the file1 / cochqtest / 9.1.2 / 364
2 / jttest / 9.2.3 / 394
3 / pagetest / 9.2.4 / 400
4 / andersontest / 9.2.4 / 402
5 / durbin2 / 9.2.4 / 404
6 / rctabm / 9.3.4 / 432
7 / rctabm2 / 9.3.4 / 433
8 / symtest / 9.3.5 / 435
9 / smtest / 9.3.5 / 436
10 / mhortest / 9.3.5 / 437
11 / kappaest / 9.3.5 / 438
12 / moditr / 10.1.1 / 464
13 / bhmod / 10.1.1 / 466
14 / gridgeman / 10.1.1 / 470
15 / simforce0 / 10.2.2 / 498
16 / paireq2 / 10.2.3 / 500
17 / simana / 10.2.4 / 503
18 / ahsimtest2 / 10.2.5 / 504
19 / schtest2 / 10.2.5 / 505
20 / simmcnemar5 / 10.2.6 / 510
21 / bbmaxg / 10.3.2 / 520
22 / bbmoest2 / 10.3.2 / 520
23 / bbdis / 10.3.2 / 521
24 / bbdist / 10.3.2 / 522
25 / forcebbm2 / 10.3.4 / 536
26 / cbbmaxg / 10.3.4 / 536
27 / cbval / 10.3.6 / 542
28 / dmff / 10.3.6 / 544
30 / fitdm2 / 10.3.6 / 545
31 / dmtest2 / 10.3.7 / 552
32 / ferrisdatf / Exercise 10.17 / 557
33 / sfmodel / 11.2.1 / 579
34 / bradley / 11.2.2 / 592
35 / bradleyv / 11.2.2 / 592
36 / tmmodel / 11.2.3 / 604
37 / ratdel / 11.3.3 / 614
38 / MF / 11.3.3 / 614
39 / anadvn / 11.3.4 / 620
40 / sddvn / 11.3.4 / 620
41 / dstest / 11.3.4 / 624
42 / riest / 11.4.2 / 627
43 / rdtest / 11.4.3 / 628
44 / rsddr / 11.43 / 635
> cochqtest
function(x)
{
n <- dim(x)[1]
m <- dim(x)[2]
tt <- sum(x)
s <- sumatrix(x)[1:n, m + 1]
t <- sumatrix(x)[n + 1, 1:m]
qch1 <- (m - 1) * (m * sum(t^2) - tt^2)
qch2 <- m * tt - sum(s^2)
qch <- qch1/qch2
pv <- 1 - pchisq(qch, m - 1)
res <- round(c(qch, pv), 3)
res
}
> cochqtest(cochqdat4)
[1] 8.234 0.041
> cochqdat4
numeric matrix: 100 rows, 4 columns.
[,1] [,2] [,3] [,4]
[1,] 0 1 1 0
[2,] 0 1 1 0
[3,] 1 1 0 1
[4,] 1 1 1 1
[5,] 1 1 0 1
[6,] 1 1 1 1
[7,] 1 1 1 1
[8,] 1 1 1 0
[9,] 0 0 1 1
[10,] 1 1 1 1
[11,] 0 1 1 1
[12,] 1 1 1 1
[13,] 1 1 1 1
[14,] 1 1 1 1
[15,] 0 1 1 1
[16,] 0 1 1 1
[17,] 1 1 1 1
[18,] 1 1 0 0
[19,] 1 1 1 1
[20,] 0 0 0 0
[21,] 1 0 0 0
[22,] 1 0 0 1
[23,] 1 1 1 0
[24,] 1 1 1 1
[25,] 0 1 1 1
[26,] 1 1 1 1
[27,] 0 1 0 1
[28,] 0 1 1 1
[29,] 1 1 1 1
[30,] 0 0 1 1
[31,] 0 1 1 1
[32,] 1 0 1 0
[33,] 1 1 1 0
[34,] 1 0 1 1
[35,] 1 0 0 1
[36,] 1 0 1 1
[37,] 1 0 0 1
[38,] 1 1 1 0
[39,] 1 1 1 1
[40,] 0 0 1 1
[41,] 1 1 1 1
[42,] 1 1 1 1
[43,] 0 1 1 1
[44,] 0 1 1 1
[45,] 1 1 1 1
[46,] 1 1 1 1
[47,] 1 0 0 1
[48,] 0 1 1 1
[49,] 1 1 1 1
[50,] 0 0 0 0
[51,] 1 0 1 1
[52,] 1 1 1 1
[53,] 1 1 1 1
[,1] [,2] [,3] [,4]
[54,] 1 1 1 1
[55,] 1 0 1 1
[56,] 1 1 1 1
[57,] 1 1 1 1
[58,] 1 1 1 1
[59,] 1 0 1 1
[60,] 1 0 1 1
[61,] 1 1 0 1
[62,] 1 1 0 1
[63,] 1 1 1 1
[64,] 1 0 0 1
[65,] 1 0 0 1
[66,] 1 1 1 1
[67,] 1 1 1 1
[68,] 1 0 0 1
[69,] 1 0 1 0
[70,] 0 0 0 0
[71,] 1 1 0 1
[72,] 0 1 1 1
[73,] 1 0 1 1
[74,] 1 0 1 1
[75,] 1 0 1 1
[76,] 0 1 1 1
[77,] 1 1 1 1
[78,] 0 0 1 1
[79,] 1 1 1 1
[80,] 0 0 0 0
[81,] 1 0 1 0
[82,] 1 0 1 1
[83,] 1 1 1 1
[84,] 1 1 1 1
[85,] 1 1 1 0
[86,] 1 0 1 1
[87,] 0 1 1 1
[88,] 1 0 1 0
[89,] 0 0 1 1
[90,] 0 0 0 0
[91,] 1 1 0 1
[92,] 1 1 1 1
[93,] 1 1 1 1
[94,] 0 0 0 0
[95,] 1 1 1 1
[96,] 1 1 1 1
[97,] 0 0 0 0
[98,] 1 1 1 1
[99,] 0 1 1 1
[100,] 1 1 1 1
> jttest
function(x, rcol, pcol)
{
#assume x<y, u2
#Jonckheere-Terpstra test, see Hollander 1999,p202.
#x matrix,rcol is columns of ratings, pcol is column of prod.
x <- as.matrix(x)
x <- x[x[, rcol] != "NA", ]
p <- leng(x[, pcol])
pk <- length(p)
jt <- 0
for(i in 1:(pk - 1))
for(j in (i + 1):pk) {
jt <- jt + mwest(x[x[, pcol] == p[i], rcol], x[x[, pcol] == p[j], rcol])[1]
}
jt
#################
nn <- length(x[, 1])
np <- seq(1, pk)
for(i in 1:pk) {
np[i] <- length(x[x[, pcol] == p[i], 1])
}
ejt <- (nn^2 - sum(np^2))/4
ejt
#################
tg <- as.vector(table(x[, rcol]))
g <- length(tg)
tt1 <- sum(tg * (tg - 1) * (2 * tg + 5))
tt2 <- sum(tg * (tg - 1) * (tg - 2))
tt3 <- sum(tg * (tg - 1))
vj1 <- nn * (nn - 1) * (2 * nn + 5) - sum(np * (np - 1) * (2 * np + 5)) - tt1
vj1 <- vj1/72
vj2 <- sum(np * (np - 1) * (np - 2))
vj2 <- vj2/(36 * nn * (nn - 1) * (nn - 2))
vj2 <- vj2 * tt2
vj3 <- sum(np * (np - 1)) * tt3
vj3 <- vj3/(8 * nn * (nn - 1))
vj <- vj1 + vj2 + vj3
vj
##############
zj <- (jt - ejt)/sqrt(vj)
zj
pv <- 1 - pnorm(abs(zj))
all <- c(jt, ejt, vj, zj, pv)
all
}
> jttest(krusdat,2,1)
[1] 16.00000000 37.50000000 88.89194139 -2.28038021 0.01129257
> krusdat
[,1] [,2]
[1,] 1 2.58
[2,] 1 2.57
[3,] 1 2.58
[4,] 1 2.76
[5,] 1 2.32
[6,] 2 2.46
[7,] 2 2.58
[8,] 2 2.49
[9,] 2 2.87
[10,] 2 2.56
[11,] 3 2.08
[12,] 3 2.30
[13,] 3 2.21
[14,] 3 2.47
[15,] 3 2.50
> pagetest
function(x)
{
#page (1963) JASA 58, 216-230
#x is matrix with n rows
#and k columns (k treatments)
n <- dim(x)[1]
k <- dim(x)[2]
rr <- sumatrix(x)[n + 1, 1:k]
pg <- sum(seq(1, k) * rr)
epg <- (n * k * (k + 1)^2)/4
vpg <- (n * (k - 1) * k^2 * (k + 1)^2)/144
z <- (pg - epg)/sqrt(vpg)
pv <- 1 - pnorm(abs(z))
all <- c(z, pv)
all <- round(all, 3)
all
res <- c(pg, epg, vpg)
res
all
}
> pagetest(pagedat2)
[1] 2.200 0.014
> pagedat2
A B C
1 2 3 1
2 2 3 1
3 2 3 1
4 1 3 2
5 1 3 2
6 3 1 2
7 1 3 2
8 3 1 2
9 1 3 2
10 1 3 2
11 1 3 2
12 1 3 2
13 1 3 2
14 1 3 2
15 3 1 2
16 1 3 2
17 1 3 2
18 1 2 3
19 1 2 3
20 2 1 3
21 1 2 3
22 2 1 3
23 2 1 3
24 1 2 3
25 2 1 3
26 1 2 3
27 1 2 3
28 2 1 3
29 2 1 3
30 1 2 3
31 1 2 3
32 1 2 3
33 1 2 3
34 1 2 3
35 1 2 3
36 2 1 3
37 2 1 3
38 2 1 3
39 1 2 3
40 3 2 1
41 3 1 2
42 3 2 1
43 3 2 1
44 3 2 1
45 3 1 2
46 3 2 1
47 3 1 2
48 3 2 1
49 3 2 1
50 3 1 2
> andersontest
function(x)
{
#This is for extended analysis for ranked data, see Best(1993).
a <- length(x[, 1])
bb <- 2 * sqrt(3/(a^2 - 1))
nn <- sum(x[1, ])
cc <- 6 * sqrt(5/((a^2 - 1) * (a^2 - 4)))
g1 <- seq(1, a)
for(j in 1:a) {
g1[j] <- bb * ((j - 1) - (a - 1)/2)
g1
}
g2 <- seq(1, a)
for(j in 1:a) {
g2[j] <- cc * ((j - 1)^2 - (a - 1) * (j - 1) + ((a - 1) * (a - 2))/6)
g2
}
y <- matrix(0, a, 2)
for(i in 1:a) {
aa1 <- sum(x[i, ] * g1)
aa2 <- sum(x[i, ] * g2)
y[i, 1] <- sqrt((a - 1)/(nn * a)) * aa1
y[i, 2] <- sqrt((a - 1)/(nn * a)) * aa2
}
dimnames(y) <- list(c(1:a), c("Location", "Spread"))
ss1 <- round(sum(y[, 1]^2), 3)
ss2 <- round(sum(y[, 2]^2), 3)
p1 <- round(1 - pchisq(ss1, a - 1), 4)
p2 <- round(1 - pchisq(ss2, a - 1), 4)
p3 <- round(1 - pchisq(ss1 + ss2, 2 * (a - 1)), 4)
cat("\n", "Summary:", "\n", "Location:", "df=", a - 1, ",", "SS=", ss1, ",", "p-value=", p1, "\n", "Spread:", "df=", a - 1, ",", "SS=", ss2, ",", "p-value=", p2, "\n", "Total:", "df=",
2 * (a - 1), ",", "SS=", ss1 + ss2, ",", "p-value=", p3, "\n")
y
}
> andersontest(andersondat)
Summary:
Location: df= 2 , SS= 24.797 , p-value= 0
Spread: df= 2 , SS= 28.244 , p-value= 0
Total: df= 4 , SS= 53.041 , p-value= 0
Location Spread
1 -2.254174 -3.5919965
2 4.057513 3.9043440
3 -1.803339 -0.3123475
> andersondat
[,1] [,2] [,3]
[1,] 42 64 17
[2,] 31 16 76
[3,] 50 43 30
> durbin2
function(t, r, k, rr)
{
#incomplete block rank test
#Hollander&Wolf 1999, P310,317,qv:see A17
#k:t;s:k;p:r
dd1 <- (12 * (t - 1))/(r * t * (k^2 - 1))
aa <- (r * (k + 1))/2
dd2 <- sum((rr - aa)^2)
dd <- dd1 * dd2
dd
pv <- 1 - pchisq(dd, t - 1)
res <- c(dd, pv)
res
}
> durbin2(7,3,3,c(3,4,7,8,9,4,7))
[1] 13.71428571 0.03299579
> rctabm2
function(x)
{
#Goodman and Kruskal's lumbda
n <- sum(x)
r <- dim(x)[1]
cc <- dim(x)[2]
rp <- seq(r)
cp <- seq(cc)
for(i in 1:r) {
rp[i] <- sum(x[i, ])
}
for(j in 1:cc) {
cp[j] <- sum(x[, j])
}
bb <- seq(1, r)
aa <- seq(1, cc)
for(i in 1:r) {
bb[i] <- max(x[i, ])
}
for(j in 1:cc) {
aa[j] <- max(x[, j])
}
lumb <- (sum(bb) - max(cp))/(n - max(cp))
luma <- (sum(aa) - max(rp))/(n - max(rp))
lum <- (sum(aa) + sum(bb) - max(cp) - max(rp))/(2 * n - max(cp) - max(rp))
################
#variance
b0 <- cbind(cp, seq(cc))
b0 <- b0[sort.list(b0[, 1]), ]
b1 <- b0[cc, 2]
bs0 <- x[, b1]
bs1 <- 0
for(i in 1:r) {
if(bs0[i] == max(x[i, ])) {
bs1 <- bs1 + max(x[i, ])
}
else {
bs1 <- bs1 + 0
}
}
a0 <- cbind(rp, seq(r))
a0 <- a0[sort.list(a0[, 1]), ]
a1 <- a0[r, 2]
as0 <- x[a1, ]
as1 <- 0
for(j in 1:cc) {
if(as0[j] == max(x[, j])) {
as1 <- as1 + max(x[, j])
}
else {
as1 <- as1 + 0
}
}
abs1 <- 0
for(j in 1:cc)
for(i in 1:r) {
if(x[i, j] == max(x[, j]) & x[i, j] == max(x[i, ])) {
abs1 <- abs1 + x[i, j]
}
else {
abs1 <- abs1 + 0
}
}
vlumb <- ((n - sum(bb)) * (sum(bb) + max(cp) - 2 * bs1))/((n - max(cp))^3)
vluma <- ((n - sum(aa)) * (sum(aa) + max(rp) - 2 * as1))/((n - max(rp))^3)
u1 <- (max(rp) + max(cp))/n
u2 <- (sum(aa) + sum(bb))/n
u3 <- (bs1 + as1 + max(bs0) + max(as0))/n
u <- (2 - u1) * (2 - u2) * (u1 + u2 + 4 - 2 * u3) - 2 * (2 - u1)^2 * (1 - abs1/n) - 2 * (2 - u2)^2 * (1 - x[a1, b1]/n)
vlum <- u/(n * (2 - u1)^4)
all <- c(lumb, luma, lum, vlumb, vluma, vlum)
all
pp <- c(u1, u2, u3, u, abs1, x[a1, b1])
pp
all
}
> rctabm2(rctabm2dat6)
[1] 0.138461538 0.206349206 0.171875000 0.005505690 0.010198083 0.005027682
> rctabm2dat6
[,1] [,2] [,3] [,4]
[1,] 8 11 8 9
[2,] 7 15 4 1
[3,] 2 9 18 8
> symtest
function(x)
{
#Testing for symmettry. Everitt 1992, p143.
k <- length(x[, 1])
xx <- x
for(i in 1:k)
for(j in 1:k) {
xx[i, j] <- (x[i, j] + x[j, i])/2
}
ch0 <- sum((x - xx)^2/xx)
ch <- 0
for(j in 2:k)
for(i in 1:(j - 1)) {
if(x[i, j] + x[j, i] == 0) {
ch <- ch
}
else {
ch <- ch + (x[i, j] - x[j, i])^2/(x[i, j] + x[j, i])
}
}
g2 <- 0
for(j in 2:k)
for(i in 1:(j - 1)) {
if(x[i, j] + x[j, i] == 0) {
g2 <- g2
}
else {
g2 <- g2 + 2 * x[i, j] * log((2 * x[i, j])/(x[i, j] + x[j, i]))
}
}
g2 <- g2
for(i in 2:k)
for(j in 1:(i - 1)) {
if(x[i, j] + x[j, i] == 0) {
g2 <- g2
}
if(x[i, j] == 0) {
g2 <- g2
}
else {
g2 <- g2 + 2 * x[i, j] * log((2 * x[i, j])/(x[i, j] + x[j, i]))
}
}
pv <- 1 - pchisq(ch, (k * (k - 1))/2)
pv2 <- 1 - pchisq(g2, (k * (k - 1))/2)
res <- round(c(ch, pv, g2, pv2), 3)
res
}
> symtest(smdat5)
[1] 12.035 0.007 12.752 0.005
> smdat5
[,1] [,2] [,3]
[1,] 65 26 4
[2,] 7 56 12
[3,] 2 9 19
> smtest
function(x)
{
#Stuart, 1955, Biometrika 42, 412.
k <- length(x[, 1])
y <- x
rs <- sumatrix(x)[k + 1, 1:k]
cs <- sumatrix(x)[1:k, k + 1]
d <- cs - rs
for(j in 1:k)
for(i in 1:k) {
if(i == j) {
y[i, j] <- cs[i] + rs[j] - 2 * x[i, j]
}
if(i != j) {
y[i, j] <- - x[i, j] - x[j, i]
}
}