Nothing
#========================================================================
catdap1 <- function(cdata, response.names = NULL, plot = 1, gray.shade = FALSE,
ask = TRUE) {
#========================================================================
if (plot != 1 && plot != 2)
plot <- 0
# cdata # original data
n <- dim(cdata)[2] # total number of variable
nsamp <- dim(cdata)[1] # sample size
# title # variable names
title <- names(cdata)
if (is.null(title)) {
title <- rep(1:n)
title <- as.character(title)
}
catdata <- t(cdata)
for (i in 1:n)
if (is.numeric(catdata[i, ]) == FALSE && is.factor(catdata[i, ]) == FALSE)
catdata[i, ] <- as.factor(catdata[i, ])
dname <- list()
for (i in 1:n)
if (is.numeric(cdata[, i]) == FALSE)
dname[[i]] <- levels(ordered(cdata[, i]))
catdap1.out <- catdap01(catdata, title, dname, response.names)
if (plot != 0)
plot.catdap1(catdap1.out, plot, gray.shade, ask)
class(catdap1.out) <- "catdap1"
return(catdap1.out)
}
#==========================================================================
catdap1c <- function(ctable, response.names = NULL, plot = 1, gray.shade = FALSE,
ask = TRUE) {
#==========================================================================
if (plot != 1 && plot != 2)
plot <- 0
# ctable # original contingency table data
id <- dim(ctable)
n <- length(id)
y <- as.array(ctable)
yn <- length(y)
nsamp <- 0
for (i in 1:yn)
nsamp <- nsamp + y[[i]]
cdata <- array(0, dim=c(n, nsamp))
nn <- yn
x <- c(0:(yn-1))
for (i in n:1) {
nn <- nn / id[i]
nd <- 0
for (j in 1:yn) {
m <- as.integer(x[j] / nn)
if (y[[j]] != 0) for (k in 1:y[[j]]) cdata[i, nd+k] <- m + 1
nd <- nd + y[[j]]
x[j] <- x[j] %% nn
}
}
# title # variable names
title <- names(dimnames(ctable))
if (is.null(title)) {
title <- rep(1:n)
title <- as.character(title)
}
dname <- dimnames(ctable)
if (is.null(dname)) {
dname <- rep(list(), n)
for (i in 1:n)
dname[[i]] <- as.character(1:id[i])
}
catdap1c.out <- catdap01(cdata, title, dname, response.names)
if (plot != 0)
plot.catdap1(catdap1c.out, plot, gray.shade, ask)
class(catdap1c.out) <- "catdap1"
return(catdap1c.out)
}
#==========================================================
catdap01 <- function(cdata, title, dname, response.names) {
#==========================================================
n <- dim(cdata)[1] # total number of variable
nsamp <- dim(cdata)[2] # sample size
# response variable numbers
ier <- 0
if (is.null(response.names) ) {
ires <- c(1:n)
l <- n
} else {
l <- length(response.names)
ires <- NULL
m <- 0
for (i in 1:l)
for (j in 1:n)
if (response.names[i] == title[j]) {
m <- m + 1
ires <- c(ires, j)
}
if (m == 0)
stop(" Response variable name is wrong.")
if (m < l)
ier <- 3001
l <- m
}
# conv # variable number and recode numbers
recode <- 0
iconv <- array(0, dim=c(20, recode))
# ex # variable numbers of explanatory variables
n1 <- 0
iex <- 0
# minmax
item1 <- rep(0, n)
item2 <- rep(0, n)
for (i in 1:n) {
item1[i] <- as.integer(min(cdata[i, ]))
item2[i] <- as.integer(max(cdata[i, ]))
}
nc <- rep(0, n)
for (i in 1:n)
nc[i] <- item2[i] - item1[i] + 1
n4 <- max(nc)
if (length(dname) == 0) {
for (i in 1:n)
if (is.numeric(cdata[i, ])) {
j1 <- item1[i]
j2 <- item2[i]
dname[[i]] <- c(j1:j2)
}
} else {
for (i in 1:n)
if (is.null(dname[[i]]) == TRUE) {
j1 <- item1[i]
j2 <- item2[i]
dname[[i]] <- c(j1:j2)
}
}
# skip # variable number and code number to delete
iskip1 <- 0
iskip <- 0
z <- .Fortran(C_catdap1,
as.integer(nsamp),
as.integer(n),
as.integer(l),
as.integer(recode),
as.integer(n1),
as.integer(iskip1),
as.integer(n4),
as.integer(item1),
as.integer(item2),
as.integer(iconv),
as.integer(ires),
as.integer(iex),
as.integer(iskip),
as.integer(cdata),
nc = integer(n),
ia = integer(n*l*n4*n4),
p = double(n*l*n4*n4),
total = integer(n*n4),
aic = double(n*l),
ord = integer((n-1)*l))
nc <- z$nc
ia <- array(z$ia, dim=c(n4, n4, l, n))
p <- array(z$p, dim=c(n4, n4, l, n))
total <- array(z$total, dim=c(n4, n))
aic <- array(z$aic, dim=c(l, n))
for (i in 1:l)
aic[i, ires[i]] <- NA
ord <- array(z$ord, dim=c(l, n-1))
ia <- aperm(ia, c(2, 1, 3, 4))
p <- aperm(p, c(2, 1, 3, 4))
cross.table <- list()
for (i in 1:l) {
i1 <- 1
i2 <- nc[ires[i]]
cross.table[[i]] <- list()
for (j in 1:n) {
cross.table[[i]][[j]] <- list()
if (j != ires[i]){
cross.table[[i]][[j]]$res <- title[ires[i]]
cross.table[[i]][[j]]$exv <- title[j]
j1 <- 1
j2 <- nc[j]
cross.table[[i]][[j]]$n <- array(ia[j1:j2, i1:i2, i, j], dim=c(j2,i2),
dimnames = list(dname[[j]], dname[[ires[i]]]))
cross.table[[i]][[j]]$p <- array(p[j1:j2, i1:i2, i, j], dim=c(j2,i2))
} else {
cross.table[[i]][[j]]$res <- title[ires[i]]
cross.table[[i]][[j]]$exv <- NULL
}
} # for (j in 1:n) end
} # for (i in 1:l) end
etotal <- matrix(list(), n)
for (j in 1:n) {
j1 <- 1
j2 <- nc[j]
etotal[j] <- list(total[j1:j2, j])
}
catdap01.out <- list(tway.table = cross.table, title = title, total = etotal,
aic = aic, aic.order = ord, ier = ier)
return(catdap01.out)
}
#=================================
print.catdap1 <- function(x,...) {
#=================================
l <- dim(x$aic)[1]
n <- dim(x$aic)[2]
ires <- rep(0, l)
dl22 <- "----------------------"
dl10 <- "----------"
bl14 <- " "
bl8 <- " "
#--------------------------------------------------------------#
cat("\n < Summary of AIC's for the two-way tables >\n\n")
#--------------------------------------------------------------#
cat(" Explanatory\n")
cat(" variables")
for (i in 1:l)
cat(sprintf(" %20s", x$tway.table[[i]][[1]]$res))
cat("\n", rep(dl10, 2), rep(dl22, l), "\n", sep="")
for (jj in 1:n) {
cat(sprintf("%20s", x$title[jj]))
for (j in 1:l) {
if (is.na(x$aic[j,jj]) == TRUE) {
cat(bl14, " - ")
ires[j] <- jj
} else {
if (is.null(x$aic[j,jj]) == FALSE)
cat(bl14, sprintf("%8.2f", x$aic[j,jj]), sep = "")
}
}
cat("\n")
}
nc <- n-1
if (nc < 6) {
nc1 <- nc
nr <- 1
} else if (nc < 11) {
nc1 <- as.integer((nc+1) / 2)
nr <- 2
} else {
nc1 <- 5
nr <- 2
}
dname <- list()
item <- rep(0, n)
for (i in 1:n)
if (i != ires[1]) {
dname[[i]] <- dimnames(x$tway.table[[1]][[i]]$n)[[1]]
dname[[ires[1]]] <- dimnames(x$tway.table[[1]][[i]]$n)[[2]]
}
for (i in 1:n)
item[i] <- length(dname[[i]])
for (i in 1:l) {
i1 <- 1
i3 <- ires[i]
i2 <- item[i3]
#------------------------------------------------------------------------------#
cat("\n\n < List of explanatory variables arranged in ascending order of AIC >\n")
#------------------------------------------------------------------------------#
cat(sprintf("\n Response variable : %s\n\n", x$title[i3]))
cat(" No. Explanatory Number of A I C Difference\n")
cat(bl14, "variables categories", bl8, bl8, "of AIC\n", sep = "")
cat(rep(dl10,7), sep = "")
cat("\n")
as <- 0.0
for (j in 1:nc) {
nv <- x$aic.order[i, j]
aic <- x$aic[i, nv]
if (j != 1)
as <- aic - as
cat(sprintf("%5i%20s %8i %8.2f %8.2f\n", j,x$title[nv], item[nv],
aic, as))
as <- aic
} # for (j in 1:nc) end
#--------------------------------------------------------------------------#
cat("\n\n < Two-way tables arranged in ascending order of AIC >\n\n")
#--------------------------------------------------------------------------#
cat(sprintf("\t\t\t\t\t( %s )\n", x$title[i3]))
nsamp <- 0
ptt <- 0.0
pt <- rep(0, i2)
cat(bl8)
for (ii in i1:i2) {
cat(sprintf("%18s", dname[[i3]][ii]))
nsamp <- nsamp + x$total[[i3]][ii]
}
for (ii in i1:i2) {
pt[ii] <- x$total[[i3]][ii] / nsamp * 100.0
ptt <- ptt + pt[ii]
}
cat(" Total\n")
for (j in 1:nc) {
nv <- x$aic.order[i, j]
j1 <- 1
j2 <- item[nv]
cat(sprintf("( %s )\n", x$title[nv]))
for( k in j1:j2) {
cat(sprintf("%10s", dname[[nv]][k]))
tc <- 0
for (ii in i1:i2)
cat(sprintf("% 8i ( %5.1f )", x$tway.table[[i]][[nv]]$n[k,ii],
x$tway.table[[i]][[nv]]$p[k,ii]))
cat(sprintf("% 8i ( %5.1f )\n", x$total[[nv]][k], ptt))
}
cat(" Total")
for (ii in i1:i2)
cat(sprintf("%8i ( %5.1f )", x$total[[i3]][ii], pt[ii]))
cat(sprintf("%8i ( %5.1f )\n\n", nsamp, ptt))
} # for (j in 1:nc) end
} # for (i in 1:l) end
if (x$ier == 3001)
warning(" some response variable name is wrong")
}
#=================================================
plot.catdap1 <- function(x, plot.type, gray.shade, ask, ...) {
#=================================================
l <- dim(x$aic)[1]
n <- dim(x$aic)[2]
ires <- rep(0, l)
for (jj in 1:n)
for (j in 1:l)
if (is.na(x$aic[j,jj]) == TRUE)
ires[j] <- jj
dname <- list()
item <- rep(0, n)
for (i in 1:n)
if (i != ires[1]) {
dname[[i]] <- dimnames(x$tway.table[[1]][[i]]$n)[[1]]
dname[[ires[1]]] <- dimnames(x$tway.table[[1]][[i]]$n)[[2]]
}
for (i in 1:n)
item[i] <- length(dname[[i]])
old.par <- par(no.readonly=TRUE)
par.ask <- FALSE
for (i in 1:l) {
if (plot.type == 1)
plot.single1(dname, ires[i], item, x$aic.order[i, ], x$aic[i, ],
x$tway.table[[i]], x$total, gray.shade, old.par, par.ask)
if (plot.type == 2)
plot.single2(dname, ires[i], item, x$aic.order[i, ], x$aic[i, ],
x$tway.table[[i]], x$total, gray.shade, old.par, par.ask)
par.ask <- ask
}
i1 <- 1
i3 <- ires[1]
i2 <- item[i3]
nsamp <- 0
for (ii in i1:i2)
nsamp <- nsamp + x$total[[i3]][ii]
plot.shade(x$title, x$aic, ires, nsamp, gray.shade, old.par, ask)
}
#=================================================
plot.shade <- function(title, aic, ires, nsamp, gray.shade, old.par, ask) {
#=================================================
new.mai <- old.par$mai
new.mai[1] <- new.mai[1] * 1.5
new.mai[2] <- new.mai[2] * 1.5
new.mai[3] <- new.mai[3] * 2.0
new.mai[4] <- new.mai[4] * 1.0
response <- NULL
explanatory <- NULL
nres <- length(ires)
for (i in nres:1)
response <- c(response, title[ires[i]])
nex <- length(title)
for (i in 1:nex)
explanatory <- c(explanatory, title[i])
d1 <- dim(aic)[1]
d2 <- dim(aic)[2]
aics <- aic / nsamp
aicsmin <- 0.0
aicsmax <- 0.0
for (i1 in 1:d1)
for (i2 in 1:d2)
if (is.na(aics[i1, i2]) == FALSE) {
aicsmin <- min(aicsmin, aics[i1, i2])
aicsmax <- max(aicsmax, aics[i1, i2])
}
if (aicsmin > -0.1)
aicsmin <- -0.15
if (aicsmax < 0)
aicsmax <- 0.5
aics1 <- array(0, c(d1, d2))
for (i in d1:1)
aics1[i, ] <- aics[d1-i+1, ]
aics1 <- aperm(aics1, c(2, 1))
aics2 <- array(0, c(d2+1, d1))
aics2[1:d2, ] <- aics1[1:d2, ]
x <- (1:nrow(aics2))
y <- (1:ncol(aics2))
bp <- c(aicsmin, -0.1, -0.05, -0.01, 0, aicsmax) - 1.0e-10
mtitle1 <- "The shaded style display of all the AIC's\n\nColor level : "
mtitle2 <- "N=Number of observations\n 5 \ : \ \ \ \ \ \ \ \ \ \ \ \ \ "
mtitle3 <- "AIC/N < -0.10\ \n 4 \ : -0.10 =< AIC/N < -0.05\ \n 3 \ : -0.05 =<"
mtitle4 <- " AIC/N < -0.01\ \n 2 \ : -0.01 =< AIC/N < \ 0.0\ \ \n 1 \ :\ \ "
mtitle5 <- "0.0\ \ \ =< AIC/N\ \ \ \ \ \ \ \ \ \ \ \ \n"
main.title <- paste(mtitle1, mtitle2, mtitle3, mtitle4, mtitle5, sep="")
if (ask == TRUE) {
par(ask = TRUE, mai = new.mai, cex.axis = 0.8)
} else {
dev.new()
par(mai = new.mai, cex.axis = 0.8)
}
pcol <- brewer.pal(11, "RdGy")
if (gray.shade == TRUE) {
scol <- pcol[c(11,10,9,8,6)]
} else {
scol <- pcol[2:6]
}
image(x, y, aics2, col = scol, breaks = bp, axes = FALSE, xlab = "", ylab = "",
main = main.title, sub = "explanatory variable", cex.main = 0.8,
cex.sub = 0.8)
legend("topright", c("5","4","3","2","1"), pch = 15, col = scol, bg = "snow2",
cex = 0.8, pt.cex = 1.1)
axis(1, 1:d2, explanatory, las = 3)
axis(2, 1:d1, response, las = 1)
box()
par(old.par)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.