.packageName <- "catspec"
# ctab: oneway, twoway, multiway percentage tables
# first argument must consist of one or more factors
# or a table object (class table, xtabs, or ftable)
# digits: number of digits after the decimal (default 2)
# type: "n" for counts, "row", "column" or "total"
# for percentages (default "n")
# row.vars:
# col.vars: same usage as ftable, ignored for one- and
# two-way tables
# percentages: FALSE==> proportions are presented rather
# than percentages (default TRUE)

# comments to John Hendrickx <John_Hendrickx@yahoo.com>

ctab<-function(...,digits=2,
        type=c("n", "row", "column", "total"),
        row.vars=NULL, col.vars=NULL,
        percentages=TRUE) {
    if (attributes(...)$class=="factor") {
        # create a table if the arguments are factors
        tbl<-table(...)
    }
    else if ("table" %in% class(...) || class(...)=="ftable") {
        # the argument is a table object (table, xtabs, ftable)
        tbl<-eval(...)
    }
    else {
        stop("first argument must be either factors or a table object")
    }

    type<-match.arg(type)

    # one dimensional table,restrict choices to "n" and "total"
    if (length(dim(tbl))==1) {
        type<-ifelse(type=="n","n","total")
    }

    # if the object is an ftable, use the row.vars and col.vars
    # use numeric indices to avoid finding the omitted
    # the object must be converted to a table to get the dimensions right
    if (class(tbl)=="ftable") {
        nrowvar<-length(names(attr(tbl,"row.vars")))
        row.vars<-1:nrowvar
        col.vars<-(1:length(names(attr(tbl,"col.vars"))))+nrowvar
        tbl<-as.table(tbl)
    }

    # marginals to exclude assuming first factor is the row vaariable,
    # second factor is the column variable
    # is overridden by row.vars or col.vars
    mrg2drop<-0
    if (type=="column") {mrg2drop<-1}
    if (type=="row") {mrg2drop<-2}
    if (type=="total" && length(dim(tbl)) > 1) {mrg2drop<-c(1,2)}


    # use row.vars and col.vars to determine the
    # marginals to use when calculating percentages
    # start by translating names to variable positions
    nms<-names(dimnames(tbl))
    if (!is.null(row.vars) && !is.numeric(row.vars)) {
        row.vars<-order(match(nms,row.vars),na.last=NA)
    }
    if (!is.null(col.vars) && !is.numeric(col.vars)) {
        col.vars<-order(match(nms,col.vars),na.last=NA)
    }
    # calculate the other if only one is given
    if (!is.null(row.vars) && is.null(col.vars)) {
        col.vars<-(1:length(dim(tbl)))[-row.vars]
    }
    if (!is.null(col.vars) && is.null(row.vars)) {
        row.vars<-(1:length(dim(tbl)))[-col.vars]
    }
    # now determine the margin as the last element
    if (type=="row" && !is.null(col.vars)) {
        mrg2drop<-col.vars[length(col.vars)]
    }
    if (type=="column" && !is.null(row.vars)) {
        mrg2drop<-row.vars[length(row.vars)]
    }
    # if row.vars is given, col.vars has been determined
    if (type=="total" && !is.null(row.vars)) {
        mrg2drop<-c(col.vars[length(col.vars)],row.vars[length(row.vars)])
    }

    marg<-(1:length(dim(tbl)))[(-mrg2drop)]

    # create percentages
    if (type=="n") {
        digits<-0
    }
    else {
        tbl<-prop.table(tbl,marg)
        if (percentages) {tbl<-tbl*100}
    }


    # use ftable for more than 2 dimensions
    # (ftable doesn't work for 1 dimension,
    # and table is nicer for 2 dimensions IMHO
    if (length(dim(tbl))>2) {
        if (is.null(row.vars)) {
            # let the second variable be the column variable
            row.vars<-names(dimnames(tbl))[-2]
            # reverse the order, last variables are groups, first is row variable
            row.vars<-rev(row.vars)
        }
        tbl<-ftable(tbl,row.vars=row.vars,col.vars=col.vars)
    }

    # get the names of the column variable
    if (class(tbl)=="ftable") {
        nms<-attr(tbl,"col.vars")[[1]]
    }
    else if (length(dim(tbl))==1) {
        nms<-dimnames(tbl)[[1]]
    }
    else{
        nms<-dimnames(tbl)[[2]]
    }

    # present the (percentage) table
    wd<-max(nchar(nms),nchar(as.integer(tbl))+digits+1)
    tbl<-formatC(tbl,format="f",width=wd,digits=digits)
    tbl
}
# function to restructure a data-frame as a "person-choice" file:
# "datamat" is the name of the data-frame
# "catvar" is the response factor,
# i.e. the dependent variable in a multinomial logistic model
# In the "person-choice" file, each record of "datamat" is duplicated
# "ncat" times, where "ncat" is the number of categories of "catvar")
# The variable "id" indexes respondents
# and is used as the stratifying variable in "clogit"
# The variable "newy" indexes response options for each respondent
# The variable "depvar" equals 1 for the record
# corresponding with the respondents actual choice
# and is 0 otherwise.
# "depvar" is the dependent variable in "clogit"
# Once "depvar" has been created, the variable "catvar" is redundant
# and it's contents can be replaced by "newy"
# In "clogit", the main effects of "catvar" will now correspond with the
# intercept term of a multinomial logit model, interactions of "catvar" with
# other independent variables will correspond with their effects
mclgen <- function (datamat,catvar) {
	stopifnot(is.data.frame(datamat))
	attach(datamat)
	stopifnot(is.factor(catvar))
	ncat <- nlevels(catvar)
	perschoice<-as.data.frame(rep(datamat,ncat))
	perschoice<-reshape(perschoice,direction="long",
		varying=lapply(names(datamat),rep,ncat),
		timevar="newy")
	perschoice<-perschoice[sort.list(perschoice$id),]
	dep<-parse(text=paste("perschoice$",substitute(catvar),sep=""))
	perschoice$depvar<-ifelse(as.numeric(eval(dep))==perschoice$newy,1,0)
	perschoice[[substitute(catvar)]]<-as.factor(perschoice$newy)
	perschoice[[substitute(catvar)]]<-factor(eval(dep),labels=levels(catvar))
	perschoice
}
# calculates BIC and AIC relative to a saturated loglinear model
# rather than relative to a null model
fitmacro<-function(object) {
	stopifnot(class(object)[1]=="glm",object$family$family=="poisson",object$family$link=="log")
	ncases<-sum(object$y)
	bic<-object$deviance-object$df.residual*log(ncases)
	aic<-object$deviance-object$df.residual*2
	cat("\n","\n")
	cat("deviance:            ",formatC(object$deviance,    width = 12, digits = 3, format = "f"),"\n")
	cat("df:                  ",formatC(object$df.residual, width = 12, digits = 0, format = "f"),"\n")
	cat("bic:                 ",formatC(bic,                width = 12, digits = 3, format = "f"),"\n")
	cat("aic:                 ",formatC(aic,                width = 12, digits = 3, format = "f"),"\n")
	cat("Number of parameters:",formatC(object$rank,        width = 12, digits = 0, format = "f"),"\n")
	cat("Number of cases:     ",formatC(ncases,             width = 12, digits = 0, format = "f"),"\n")
	cat("\n","\n")
}

# utility function to check if the variables are factors
# with the same number of categories
# called by functions for mobility models below
check.square <- function(rowvar,colvar,equal=TRUE) {
	stopifnot(is.factor(rowvar))
	stopifnot(is.factor(colvar))
	if (equal) {
		stopifnot(nlevels(rowvar)==nlevels(colvar))
	}
}

# Quasi-independence
mob.qi <- function(rowvar,colvar,constrained=FALSE,print.labels=FALSE) {
	check.square(rowvar,colvar)
	if (constrained) {
		qi <- ifelse(rowvar==colvar, 1, 0)
		nms<-c("diagonal")
	}
	else {
		qi <- ifelse(rowvar==colvar, rowvar, 0)
		nms<-levels(rowvar)
	}

	qi<-factor(qi)
	qi<-C(qi,contr.treatment,base=1)
	if (print.labels) {
		levels(qi)<-c("offdiag",nms)
	}
	qi
}

# symmetric interaction effects
mob.symint <- function(rowvar,colvar,print.labels=FALSE) {
	check.square(rowvar,colvar)
	# remove factor levels to avoid messy output
	if (!print.labels) {
		attr(rowvar,"levels")<-1:nlevels(rowvar)
		attr(colvar,"levels")<-1:nlevels(colvar)
	}
	mdl<-model.matrix(~rowvar*colvar)
	intrct<-mdl[,attr(mdl,"assign")==3]
	# remove factor names
	colnames(intrct)<-sub("rowvar","",colnames(intrct))
	colnames(intrct)<-sub("colvar","",colnames(intrct))
	w<-ncol(intrct)
	x<-matrix(1:w,sqrt(w),sqrt(w))
	symint<-intrct[,t(x)[lower.tri(x,diag=TRUE)]]+intrct[,x[lower.tri(x,diag=TRUE)]]
	symint
}

# equal main effects, Hope's halfway model
mob.eqmain <- function(rowvar,colvar,print.labels=FALSE) {
	check.square(rowvar,colvar)
	if (!print.labels) {
		attr(rowvar,"levels")<-1:nlevels(rowvar)
		attr(colvar,"levels")<-1:nlevels(colvar)
	}
	rmat<-model.matrix(~rowvar)
	rmat<-rmat[,attr(rmat,"assign")==1]
	cmat<-model.matrix(~colvar)
	cmat<-cmat[,attr(cmat,"assign")==1]
	eqmain<-rmat+cmat
	colnames(eqmain)<-sub("rowvar","",colnames(eqmain))
	eqmain
}
# Crossings parameter
mob.cp <- function(rowvar,colvar) {
	check.square(rowvar,colvar)
	cp<-NULL
	rvar<-as.numeric(rowvar)
	cvar<-as.numeric(colvar)
	for (i in 1:(nlevels(rowvar)-1)) {
		cp<-cbind(cp,as.numeric((rvar <= i & cvar > i) | (rvar > i & cvar <= i)))
	}
	cp
}

# Uniform association
mob.unif <- function(rowvar,colvar) {
	check.square(rowvar,colvar,equal=FALSE)
	as.numeric(rowvar)*as.numeric(colvar)
}

# RC model 1 (unequal row and column effects, page 58)
# Fits a uniform association parameter and row and column effect
# parameters. Row and column effect parameters have the
# restriction that the first and last categories are zero.
mob.rc1 <- function(rowvar,colvar,equal=FALSE,print.labels=FALSE) {
	# the number of row and column categories need not be equal
	# unless an equality restriction is applied
	check.square(rowvar,colvar,equal=equal)
	# use numbers rather than factor levels by default
	if (!print.labels) {
		attr(rowvar,"levels")<-1:nlevels(rowvar)
		attr(colvar,"levels")<-1:nlevels(colvar)
	}

	# row effects, first and last category constrained to 0
	# multiplied by column variable as continuous
	rowvar<-C(rowvar,contr.treatment,base=1)
	rmat<-model.matrix(~rowvar)
	rmat<-rmat[,attr(rmat,"assign")==1]
	rmat<-rmat[,1:(ncol(rmat)-1)]*as.numeric(colvar)

	# column effects, same construction
	colvar<-C(colvar,contr.treatment,base=1)
	cmat<-model.matrix(~colvar)
	cmat<-cmat[,attr(cmat,"assign")==1]
	cmat<-cmat[,1:(ncol(cmat)-1)]*as.numeric(rowvar)

	u<-mob.unif(rowvar,colvar)
	# add a name for the uniform association effect
	dim(u)<-c(length(u),1)
	colnames(u)<-c("U")

	if (equal) {
		# rowmain and colmain are added to impose an equality restriction
		colnames(rmat)<-sub("rowvar","RC",colnames(rmat))
		rc1<-cbind(rmat+cmat,u)
	}
	else {
		colnames(rmat)<-sub("rowvar","R",colnames(rmat))
		colnames(cmat)<-sub("colvar","C",colnames(cmat))
		rc1<-cbind(rmat,u,cmat)
	}
	rc1
}
