Nothing
# -------------------------------------------------------------------------------------------------------------------------- #
# -------------------------------------------------------------------------------------------------------------------------- #
#wijfmat <- matrix( c(
#92, 107, 247,
#255, 89, 89,
#92, 203, 92,
#255, 177, 17,
#170, 97, 187,
#255,255,95,
#255,137,235,
#145, 101, 62,
#193,193,193,
#92, 229, 214,
#201, 255, 135,
#255, 224, 204,
#173,45,92,
#227, 196, 239,
#226, 212, 149,
#204, 241, 255,
#87, 142, 82
#), ncol=3, byrow = TRUE)
#wf <- apply(wijfmat,1,function(z){
# rgb(z[1]/255,z[2]/255,z[3]/255)
#})
wijffelaars <- c("#5C6BF7", "#FF5959", "#5CCB5C", "#FFB111", "#AA61BB", "#FFFF5F", "#FF89EB", "#91653E", "#C1C1C1", "#5CE5D6", "#C9FF87",
"#FFE0CC", "#AD2D5C", "#E3C4EF", "#E2D495", "#CCF1FF", "#578E52")
draw = function(H,W,X,Y,alpha = 1,border = "black",bg = "white",vp=NULL, lwd = 1){
grid.rect(x = unit(X, "npc"), y = unit(Y, "npc"),
width = unit(W, "npc"), height = unit(H, "npc"),
just = c("left","bottom"),
default.units = "npc", name = NULL,
gp=gpar(col=border,fill = bg, alpha = alpha, lwd = lwd), draw = TRUE, vp = vp
)
}
rmb = function(x,...){
UseMethod("rmb")
}
rmb.table = function(x, col.vars = NULL, spine = FALSE, circular = FALSE, eqwidth = FALSE, cat.ord = NULL,
freq.trans = NULL, max.scale = 1, use.na = FALSE, expected = NULL, residuals = NULL,
model.opt = list(), gap.prop = 0.2, gap.mult = 1.5, col = "hcl",col.opt = list(), label = TRUE,label.opt = list(),vp = NULL,...){
#vn <- names(dimnames(x))
#if(is.null(vn)){
# vn <- LETTERS[1:length(dim(x))]
#}
#if(is.numeric(col.vars)){
# vn <- c(vn[-col.vars],vn[col.vars])
# col.vars <- sort(seq_along( vn ) %in% col.vars)
#}
#formula <- as.formula(paste("~",paste(vn,collapse="+"),sep=""))
#data <- as.data.frame(ftable(x))
nd <- length(dim(x))
data <- subtable(x,1:nd)
formula <- as.formula(paste("~",paste(names(data)[1:nd],collapse="+"),sep=""))
rmb.formula(formula, data, col.vars = col.vars, spine = spine, circular = circular,
eqwidth = eqwidth, cat.ord = cat.ord, freq.trans = freq.trans, max.scale = max.scale,
use.na = use.na, expected = expected, residuals = residuals, model.opt = model.opt,
gap.prop = gap.prop, gap.mult = gap.mult, col = col,col.opt = col.opt, label = label,
label.opt = label.opt, vp = vp)
}
rmb.ftable = function(x, col.vars = NULL, spine = FALSE, circular = FALSE, eqwidth = FALSE, cat.ord = NULL,
freq.trans = NULL, max.scale = 1, use.na = FALSE, expected = NULL, residuals = NULL,
model.opt = list(), gap.prop = 0.2, gap.mult = 1.5, col = "hcl",col.opt = list(), label = TRUE,label.opt = list(),vp = NULL,...){
rv <- names(attr(x,"row.vars"))
cv <- names(attr(x,"col.vars"))
col.vars <- NULL
#if(is.null(col.vars)){
col.vars <- c(rep(F,length(rv)),rep(T,length(cv)))
#}
formula <- as.formula(paste("~",do.call(paste,c(as.list(c(rv,cv)),sep="+")),sep=""))
data <- as.data.frame(x)
rmb.formula(formula, data, col.vars = col.vars, spine = spine, circular = circular,
eqwidth = eqwidth, cat.ord = cat.ord, freq.trans = freq.trans, max.scale = max.scale,
use.na = use.na, expected = expected, residuals = residuals, model.opt = model.opt,
gap.prop = gap.prop, gap.mult = gap.mult, col = col,col.opt = col.opt, label = label,
label.opt = label.opt,vp = vp)
}
rmb.glm = function(x, col.vars = NULL, spine = FALSE, circular = FALSE, eqwidth = FALSE, cat.ord = NULL,
freq.trans = NULL, max.scale = 1, use.na = FALSE, expected = NULL, residuals = NULL,
model.opt = list(), gap.prop = 0.2, gap.mult = 1.5, col = "hcl",col.opt = list(), label = TRUE,label.opt = list(),vp = NULL,...){
model.family <- x$family$family
stopifnot(model.family %in% c("poisson","binomial"))
if(model.family == "poisson"){
fi <- which(names(x$data) == "Freq")
tms <- terms(x$formula)
vars <- attr(tms,"term.labels")[ which(attr(tms,"order") == 1)]
formula <- paste( toString(x$formula[2]),"~", paste(vars,collapse="+"),sep="" )
if("resid.type" %in% names(model.opt)){
rtp <- model.opt$resid.type
}else{
rtp <- "pearson"
}
data <- subtable(x$data, match(vars,names(x$data)[-fi]), keep.zero=TRUE)
x <- update(x, data = data)
residuals <- residuals(x,type=rtp)
expected <- x$fitted
}
if(model.family == "binomial"){
tms <- terms(x$formula)
vars <- attr(tms,"term.labels")[ which(attr(tms,"order") == 1)]
vars <- c(vars[-1],vars[1])
formula <- paste( "Freq ~ ",toString(x$formula[2]), paste(vars,collapse="+"),sep="" )
nv <- length(vars)+1
unweighted <- all(x$prior.weights==1)
data <- subtable(x$data, 1:nv, keep.zero=TRUE, allfactor = TRUE)
if(unweighted){
x <- update(x, data = data, weights = data$Freq)
}else{
x <- update(x, data = data)
}
if("resid.type" %in% names(model.opt)){
rtp <- model.opt$resid.type
}else{
rtp <- "pearson"
}
residuals <- residuals(x,type=rtp)
expected <- x$fitted
}
rmb.formula(formula, data, col.vars = col.vars, spine = spine, circular = circular,
eqwidth = eqwidth, cat.ord = cat.ord, freq.trans = freq.trans, max.scale = max.scale,
use.na = use.na, expected = expected, residuals = residuals, model.opt = model.opt,
gap.prop = gap.prop, gap.mult = gap.mult, col = col,col.opt = col.opt, label = label,
label.opt = label.opt, vp = vp)
}
#rmb.polr = function(x, col.vars = NULL, spine = FALSE, circular = FALSE, eqwidth = FALSE, cat.ord = NULL,
# freq.trans = NULL, max.scale = 1, use.na = FALSE, expected = NULL, residuals = NULL,
# model.opt <- list(), gap.prop = 0.2, gap.mult = 1.5, col = "hcl",col.opt = list(), label = TRUE,label.opt = list(),...){
#
# model.family <- x$family$family
# stopifnot(model.family %in% c("poisson","binomial"))
#
#
# vars <- dimnames(attr(x$terms,"factors"))[[1]]
# vars <- c(vars[-1],vars[1])
# nv <- length(vars)
# formula <- paste( "Freq~", paste(vars,collapse="+"),sep="" )
#
# data <- x$model
#
# unweighted <- ncol(data) == nv
#
# if(unweighted){
# data <- subtable(data, 1:nv, keep.zero=FALSE, allfactor = TRUE)
# x <- update(x, data = data, weights = data$Freq)
# }else{
# names(data)[which(!(names(data) %in% vars))] <- "Freq"
# }
#
# residuals <- residuals(x,type=resid.type)
# expected <- x$fitted
#
# rmb.formula(formula, data, col.vars = col.vars, spine = spine, circular = circular,
# eqwidth = eqwidth, cat.ord = cat.ord, freq.trans = freq.trans, max.scale = max.scale,
# use.na = use.na, expected = expected, residuals = residuals, model.opt = model.opt,
# gap.prop = gap.prop, gap.mult = gap.mult, col = col,col.opt = col.opt, label = label,
# label.opt = label.opt)
#}
rmb.formula = function(formula, data, col.vars = NULL, spine = FALSE, circular = FALSE, eqwidth = FALSE, cat.ord = NULL,
cut = NULL, innerval = 1, freq.trans = NULL, num.mode = FALSE, max.scale = 1, use.na = FALSE, expected = NULL, residuals = NULL,
model.opt = list(), gap.prop = 0.2, gap.mult = 1.5, col = "hcl",col.opt = list(), label = TRUE,label.opt = list(), vp = NULL,...){
# use.expected.values = FALSE, mod.type = "poisson", resid.type = "pearson",
# resid.display = "both",max.rat = 2.0, resid.max = NULL, cut.rs = 5
# boxes = TRUE, lab.tv = FALSE, varnames = TRUE, abbrev = FALSE, lab.cex = 1.2, yaxis = TRUE
if(!is.null(vp)){
if(!inherits(vp,"viewport")){
stopifnot(length(vp)>1)
if(is.null(current.viewport()$layout)){
print("No layout specified. Try:")
print("mat.layout <- grid.layout(nrow, ncol); grid.newpage(); vp.mat <- viewport(layout = mat.layout);pushViewport(vp.mat)")
}
vp=viewport(layout.pos.row = vp[1], layout.pos.col = vp[2])
}
pushViewport(vp)
}else{
grid.newpage()
}
min.alpha <- 0.1
if(!is.null(expected)){
if(is.logical(expected)){
if(expected == FALSE){
expected <- NULL
}
}
}
# ----- get model parameters from model.opt or use defaults ------ #
if( "use.expected.values" %in% names(model.opt) ){
use.expected.values <- model.opt$use.expected.values
}else{
use.expected.values <- FALSE
}
if( "mod.type" %in% names(model.opt) ){
mod.type <- model.opt$mod.type
}else{
mod.type <- "poisson"
}
if( "resid.type" %in% names(model.opt) ){
resid.type <- model.opt$resid.type
}else{
resid.type <- "pearson"
}
if( "resid.display" %in% names(model.opt) ){
resid.display <- model.opt$resid.display
}else{
resid.display <- "both"
}
if( "max.rat" %in% names(model.opt) ){
max.rat <- model.opt$max.rat
}else{
max.rat <- 2.0
}
if( "resid.max" %in% names(model.opt) ){
resid.max <- model.opt$resid.max
}else{
resid.max <- NULL
}
if( "cut.rs" %in% names(model.opt) ){
cut.rs <- model.opt$cut.rs
}else{
cut.rs <- 5
}
if( "bg" %in% names(col.opt) ){
bg <- col.opt$bg
}else{
bg <- NA #hsv(0,0,0, alpha=0.02)#hcl(240,10,85)
}
if( "col2" %in% names(col.opt) ){
col2 <- col.opt$col2
}else{
col2 <- alpha("grey",0.25)
}
if( "ang" %in% names(col.opt) ){
ang <- col.opt$ang
}else{
ang <- 360
}
if( "bgs" %in% names(col.opt) ){
bgs <- col.opt$bgs
}else{
bgs <- hsv(0,0,0, alpha=0.02)
}
if( "alpha.r" %in% names(col.opt) ){
stopifnot( col.opt$alpha.r <= 1 & col.opt$alpha.r > 0 )
alpha.r <- -log(col.opt$alpha.r)
}else{
alpha.r <- -log(0.2)
}
if( "line.col" %in% names(col.opt) ){
line.col <- col.opt$line.col
}else{
if(!is.null(expected)){
if(spine | circular){
line.col <- NA
}else{
line.col <- alpha(1,0.9)
}
}else{
line.col <- ifelse( !spine & !circular, "darkgrey", NA)
}
}
if( "rect" %in% names(col.opt) ){
rect <- col.opt$rect
}else{
rect <- bgs
}
# ----- get label parameters from label.opt or use defaults ------ #
if( "yaxis" %in% names(label.opt) ){
yaxis <- label.opt$yaxis
}else{
yaxis <- TRUE & !num.mode
}
if( "boxes" %in% names(label.opt) ){
boxes <- label.opt$boxes
}else{
boxes <- TRUE
}
if( "lab.tv" %in% names(label.opt) ){
lab.tv <- label.opt$lab.tv
}else{
lab.tv <- FALSE
}
if( "varnames" %in% names(label.opt) ){
varnames <- label.opt$varnames
}else{
varnames <- TRUE
}
if( "abbrev" %in% names(label.opt) ){
abbrev <- label.opt$abbrev
}else{
abbrev <- NULL
}
if( "lab.cex" %in% names(label.opt) ){
lab.cex <- label.opt$lab.cex
}else{
lab.cex <- 1.2
}
if( "s0" %in% names(label.opt) ){
s0 <- label.opt$s0
#s1 <- label.opt$s[2]
#s2 <- label.opt$s[3]
#s4 <- label.opt$s[4]
}else{
s0 <- 0.02 # border
}
s1 <- 0.04 # yaxis
s2 <- 0.06 # labs
s3 <- 0.1 # expected scale
tmp <- which(cut == FALSE)
if(any(tmp)) cut[tmp] <- 0
tmp <- which(cut == TRUE)
if(any(tmp)) cut[tmp] <- 12
# ----- parameter check 1 ----- #
dset <- data.frame(data)
f <- as.formula(formula)
cut.rv <- TRUE
if(max.scale > 1 ){
stop(simpleError("Wrong max.scale specification! (0,1]"))
}
if(!( (gap.prop < 1) & (gap.prop > 0) ) ){
stop(simpleError("Wrong gap.prop specification! (0,1)"))
}
if( !((min.alpha <= 1) & (min.alpha >= 0)) ){
stop(simpleError("Wrong alpha specification! [0,1]"))
}
stopifnot(is.logical(label)|is.numeric(label))
#if( !("base.alpha" %in% ls()) ){
# base.alpha <- 1
#}
# ----- terms extraction ----- #
fterms <- attr(terms(f),"term.labels")
nv <- length(fterms)
# ----- cutting numeric variables ----- #
int.vars <- sapply(fterms, function(z){
v <- data[,match(z,names(data))]
if(is.numeric(v)){
all( (v%%1)==0 )
}else{
FALSE
}
})
if(length(cut) > 1){
int.vars[cut>0] <- FALSE
}
# integer variables will not be cut unless the value of cut at this particular position is positive.
# cut variable i iff cut[i] > 0 and Vi is numeric.
num.vars <- sapply(fterms, function(z) is.numeric(data[,match(z,names(data))]))
cuttf <- (cut > 0) & ( !int.vars & num.vars )
if(is.null(cut)){
cut <- rep(12,nv)
}
if(length(cut) == 1){
cut <- rep(cut, nv)
}
if(length(cut) < nv){
cut <- c(cut, rep(0,nv-length(cut)))
}
for(i in fterms[cuttf]){
j <- which(names(data) == i)
if(innerval < 1){
iv <- innerval(dset[,j],p=innerval)
too.low <- which(dset[,j] < iv[1])
too.high <- which(dset[,j] > iv[2])
#data[c(too.high,too.low),j] <- NA
dset[c(too.high,too.low),j] <- NA
}
dset[,j] <- cut(dset[,j], breaks = cut[match(i, fterms)], right = FALSE)
levels(dset[,j]) <- paste(colsplit(levels(dset[,j]),",",1:2)[,1],")", sep="")
}
#if(innerval < 1 & any(num.vars)){
#kill <- sapply(dset[,match(fterms,names(data))], is.na)
#kill <- apply(kill,1,any)
#if(any(kill)){
# dset <- dset[-which(kill),]
# dims of data and dset are now potentially different!
#}
#}
if(!use.na){
kill <- sapply(dset[,match(fterms,names(data))], is.na)
kill <- apply(kill,1,any)
#dset <- na.omit(dset)
if(any(kill)){
dset <- dset[-which(kill),]
}
}else{
for( i in 1:(ncol(dset)-1) ){
ind <- is.na(dset[,i])
if(any(ind)){
levels(dset[,i]) <- c(levels(dset[,i]),"N/A")
dset[which(ind),i] <- "N/A"
}
}
}
if(is.numeric(col.vars)){
col.vars <- seq_along( fterms ) %in% col.vars
}
#row.vars <- seq_along(fterms)[-c(col.vars)]
if( !is.null(expected) ){
both <- resid.display %in% sapply(1:4, function(z) abbreviate("both",z))
res.val.only <- resid.display %in% c(sapply(1:6, function(z) abbreviate("values",z)),sapply(1:6, function(z) abbreviate("value",z)))
disp.res <- both | res.val.only
if( disp.res & all( sapply(expected, function(z) all(z == TRUE))) ){
if( mod.type == "poisson" ){
expected <- list(1:(nv-1),nv)
}else{
expected <- as.list(1:(nv-1))
}
}
if( resid.display %in% sapply(1:6, function(z) abbreviate("values",z)) ){
col <- "rgb"
col.opt$alpha <- 0.3
}
}else{
res.val.only <- FALSE
}
nv0 <- nv
if(is.null(col.vars)){
col.vars <- rep(c(T,F),nv)[1:nv]
}
col.vars[nv] <- TRUE
tv <- fterms[nv]
stopifnot(is.logical(col.vars) & length(col.vars) == nv)
ind <- match(fterms,names(dset))
#if(num.mode){
# dset.num <- as.data.frame(data)[,ind]
#}
# >>> a dummy-variable is used to handle the case of no row-variables
if( sum(!col.vars) < 1 ){
dset <- data.frame(probability = rep("prob",nrow(dset)),dset)
nv<-nv+1
col.vars <- c(F,col.vars)
ind <- c(1,ind+1)
}
if( sum(col.vars) < 2 ){
dset <- data.frame(distribution = rep(tv,nrow(dset)),dset)
nv<-nv+1
col.vars <- c(T,col.vars)
ind <- c(1,ind+1)
if("probability" %in% names(dset)){ expected <- NULL }
}
# >>> preparing the labels for plotting including the abbreviation option
if(is.numeric(label)){
lab <- TRUE
tmp <- rep(FALSE,nv)
tmp[label] <- TRUE
label <- tmp
}else{
lab <- any(label)
if(length(label) == 1 & lab){
label <- rep(TRUE,nv)
}
if(!lab){
label <- rep(FALSE,nv)
}
}
if(length(abbrev) < length(ind)){
abbrev <- rep(abbrev,length(ind))[seq_along(ind)]
}
if(is.null(abbrev)){
rclabs <- lapply(dset[,ind],function(x) levels(as.factor(x)))
}else{
#rclabs <- lapply(dset[,ind],function(x) abbreviate(levels(as.factor(x)),abbrev))
rclabs <- mapply(function(y,z) abbreviate(levels(y),z),y = dset[,ind], z = as.list(abbrev),SIMPLIFY = FALSE)
}
lab.tv <- ifelse(spine,FALSE,lab.tv)
label[nv] <- lab.tv
if( sum(label) == 0 ){
lab <- FALSE
}
rlabs <- lapply(which( (!col.vars)*label > 0),function(x) rclabs[[x]])
clabs <- lapply(which(col.vars*label > 0),function(x) rclabs[[x]])
# orig.labs <- lapply(dset[,ind],function(x) levels(as.factor(x)))
if( length(terms(f)) > 2 ){
fi <- which( suppressWarnings(names(dset) == terms(f)[[2]]))
if("Freq" %in% names(dset)[-fi]){
names(dset)[ which(names(dset) == "Freq") ] <- "Freq2"
}
names(dset)[fi] <- "Freq"
}
dset <- subtable(dset,ind,keep.zero=FALSE,allfactor=TRUE)
# ----- predefined expected values ----- #
if( is.numeric(expected) ){
stopifnot( length(expected) == length(dset$Freq) )
mod.type <- "gen"
pred <- expected
}else{
if(suppressWarnings(max(unlist(expected))) > (nv - (mod.type=="logit")) ){
stop(simpleError("Wrong expected specification!"))
}
}
if(is.null(cat.ord)){
ntc <- nlevels(dset[,nv])
cat.ord <- 1:ntc
}
if(length(cat.ord) > 1 & !spine){
levels(dset[,nv])[which(! 1:nv %in% cat.ord )] <- NA
dset[,nv] <- factor(dset[,nv],levels=levels(dset[,nv])[rank(cat.ord)])
if(lab.tv){
clabs[[length(clabs)]] <- levels(dset[,nv])
}
}
if(!use.na){
dset <- na.omit(dset)
}else{
for( i in 1:(ncol(dset)-1) ){
ind <- is.na(dset[,i])
if(any(ind)){
levels(dset[,i]) <- c(levels(dset[,i]),"N/A")
dset[which(ind),i] <- "N/A"
}
}
}
# ----- descriptive parameters ----- #
ntc <- nlevels(dset[,nv])
ntc0 <- ntc
#nlvl <- sapply(dset[,(sapply(dset,class)=="factor")], nlevels)
nlvl <- sapply(dset[,-ncol(dset)], nlevels)
col.nlvl <- nlvl[which(col.vars)]
row.nlvl <- nlvl[which(!col.vars)]
nc <- prod(col.nlvl)
nr <- prod(row.nlvl)
ncv <- length(col.nlvl)
nrv <- length(row.nlvl)
# ----- computation of the underlying (relative) frequencies ----- #
tt1 <- ftable(tapply(dset$Freq,as.list(dset[,1:nv]),sum),col.vars=which(col.vars))
tt2 <- spread(ftable(tapply(dset$Freq,as.list(dset[,1:(nv-1)]),sum),col.vars=which(col.vars[1:(nv-1)])),ncol=ntc)
tt1[is.na(tt1)] <- 0
tt2[is.na(tt2)] <- 1
H0 <- tt1/tt2
H <- H0
# ----- a few more auxiliary variables ----- #
rind <- which( (!col.vars)*label > 0)
cind <- which(col.vars*label > 0)
nrl <- length(rind)*lab
ncl <- length(cind)*lab
num.mode.x <- num.mode
num.mode.y <- num.mode
col.gap.prop <- ifelse(nc > ntc, gap.prop * min(1,ncv/nrv), 0)
row.gap.prop <- ifelse(nr > 1, gap.prop * min(1,nrv/ncv), 0)
if( is.factor(data[, ind[cind[length(cind)]] ])){
num.mode.x <- FALSE
}
if( is.factor(data[, ind[rind[length(rind)]] ]) ){
num.mode.y <- FALSE
}
if(num.mode.x & num.mode.y){
gap.prop <- 0.001
}
if(num.mode.y){
row.gap.prop <- 0.001
}
if(num.mode.x){
col.gap.prop <- 0.001
}
# ----- model computation in expected mode ----- #
if(!is.null(expected) ){
int.dset <- lapply(dset,as.integer)
dset <- dset[do.call("order",int.dset),]
if(!(mod.type %in% c("poisson","polr", "gen"))){
stop(simpleError("Wrong mod.type specification!"))
}
# ------ logit response residuals ----- #
if(mod.type == "polr"){
if(! (resid.type %in% c("response","prob","prop","rat"))){
print(simpleWarning("Argument resid.type ignored. Only response residuals are implemented for polr models."))
resid.type<-"response"
}
nameslist <- names(dset)[1:(nv-1)]
nameslist <- nameslist[which( !(nameslist %in% c("probability","distribution")) )]
single.terms <- paste( nameslist ,collapse = "+")
interaction.terms <- do.call("paste",c(lapply(expected,function(x) do.call("paste",c(as.list(nameslist[x]),sep="*"))),sep="+"))
full.terms <- paste(nameslist,collapse = "*")
mod.formula <- as.formula(paste(fterms[nv0],"~",single.terms,"+",interaction.terms,sep=""))
modS.formula <- as.formula(paste(fterms[nv0],"~",full.terms,sep=""))
if (requireNamespace("MASS", quietly = TRUE)) {
modS <- MASS::polr( formula = as.formula(modS.formula),data = dset, weights = dset$Freq, method ="logistic")
mod <- MASS::polr( formula = as.formula(mod.formula),data = dset, weights = dset$Freq, method ="logistic")
}
pred <- predict(mod,newdata = subtable(dset,c(1:(nv-1)),keep.zero=F),type="probs")
fit.values <- as.vector(t(pred))
tt1r <- ftable(tapply(fit.values,as.list(dset[,1:nv]),sum),col.vars=which(col.vars))
resid.mat <- (H - tt1r)
p.value<-anova(mod,modS)[2,7]
if( use.expected.values ){
H0 <- tt1r#/tt2
H00 <- H
H <- H0
resid.mat <- -resid.mat
}else{
H00 <- tt1r#/tt2
# H <- H
}
}
# ----- poisson model ----- #
if(mod.type == "poisson"){
nameslist <- names(dset)[1:(nv)]
nameslist <- nameslist[which( !(nameslist %in% c("probability","distribution")) )]
single.terms <- do.call("paste",c(as.list( nameslist ),sep="+"))
interaction.terms <- do.call("paste",c(lapply(expected,function(x) do.call("paste",c(as.list(nameslist[x]),sep="*"))),sep="+"))
mod.formula <- paste("Freq~",single.terms,"+",interaction.terms,sep="")
mod <- glm( as.formula(mod.formula),data = dset, family = "poisson")
if(resid.type %in% c("prob","rat","prop")){
mod.residuals <- residuals(mod,type="response")
resid.mat <- ftable(xtabs(mod.residuals~.,data=dset[,1:nv]),col.vars=which(col.vars))
resid.mat <- resid.mat / tt2
resid.type <- "response"
}else{
mod.residuals <- residuals(mod,type=resid.type)
resid.mat <- ftable(xtabs(mod.residuals~.,data=dset[,1:nv]),col.vars=which(col.vars))
}
#resid.mat <- ftable(tapply(mod.residuals,as.list(dset[,1:nv]),sum),col.vars=which(col.vars))
Df <- mod$df.residual
Dev <- mod$deviance
p.value <- 1-pchisq(Dev,Df)
pred <- predict(mod,type="response")
#pred.mat <- ftable(tapply(pred,as.list(dset[,1:nv]),sum),col.vars=which(col.vars))
pred.mat <- ftable(xtabs(pred~.,data=dset[,1:nv]),col.vars=which(col.vars))
if( use.expected.values ){
H0 <- pred.mat/tt2
H00 <- H
H <- H0
resid.mat <- -resid.mat
}else{
H00 <- pred.mat/tt2
}
}
if(mod.type == "gen"){
p.value <- 1
# pred is already defined by expected = c(...)
pred.mat <- ftable( xtabs(pred~., data = data[,1:nv]), col.vars = which(col.vars))
if(is.null(residuals)){
mod.residuals <- dset$Freq - pred
resid.type <- "response"
}else{
stopifnot(length(residuals) == length(pred))
mod.residuals <- residuals
}
resid.mat <- ftable( xtabs(mod.residuals~., data = data[,1:nv]), col.vars = which(col.vars))
if( use.expected.values ){
H0 <- pred.mat/tt2
H00 <- H
H <- H0
resid.mat <- -resid.mat
}else{
H00 <- pred.mat/tt2
}
}
}else{
disp.res <- FALSE
H00 <- H0
}
# ----- modify parameters for rmbplot mode ----- #
if( spine | circular ){
nlvl[nv] <- 1
nc <- nc/ntc
ntc <- 1
}
# preparing the transformation parameters
if( is.list(freq.trans) ){
stopifnot(freq.trans[[1]] == "sqrt")
stopifnot(freq.trans[[2]] > 0)
tfreq <- "sqrt"
tfn <- freq.trans[[2]]
}else{
tfreq <- freq.trans
check <- FALSE
if(is.null(tfreq)){
tfreq <-"sqrt"
tfn <- 1
check <- TRUE
}else{
if( tfreq == "sqrt" ){
tfn <- 2
check <- TRUE
}
if( tfreq == "log" ){
tfn <- exp(1)
check <- TRUE
}
if( tfreq %in% c("1","c","const","eq","equal") ){
tfreq <- "const"
tfn <- 1
check <- TRUE
}
}
if(!check){
stop(simpleError("Wrong freq.trans specification!"))
}
}
# ----- computation of rectangle widths, heights and x/y coordinates in matrix form ----- #
W <- matrix(1,ncol=nc,nrow=nr)
S <- space(nlvl,col.vars,gap.prop=gap.prop,gap.mult=gap.mult,last.col.zero = TRUE,last.row.zero = FALSE)
# better usage of space here for the num.mode?
x <- ( seq(0,1-1/nc,1/nc) )*(1-col.gap.prop)
if(nc*(nlevels(dset[,nv])^(spine|circular)) > nlevels(dset[,nv])){
x <- x + c(0,cumsum(S[[1]]))*col.gap.prop/sum(S[[1]])
}
if( suppressWarnings(any(S[[2]])) ){
y <- rev((seq(0,1-1/nr,1/nr))*(1-row.gap.prop) + c(0,cumsum(S[[2]]))*row.gap.prop/sum(S[[2]])) #rev test
}else{
y <- 0
}
X <- spread(t(matrix( x )),nrow=nr)
Y0 <- spread(matrix( y ),ncol=nc)
Y <- Y0
X2 <- X[,seq(1,nc,ntc)]
Y2 <- Y[,seq(1,nc,ntc)]
H2 <- matrix(1,ncol=nc/ntc,nrow=nr)
tt3 <- ftable(tapply(dset$Freq,as.list(dset[,1:(nv-1)]),sum),col.vars=which(col.vars[1:(nv-1)]))
tt3[is.na(tt3)] <- 0
W2 <- tt3 / max(tt3)
if(nc > nlevels(dset[,nv])){
width.cor <- (ntc-1)*S[[1]][1]/sum(S[[1]])*col.gap.prop
}else{
width.cor <- 0
}
if( eqwidth ){
W2trans <- W2
if( tfreq == "const" ){
W2trans <- 1*(W2>0)
}
if( tfreq == "log" ){
W2trans <- log(tt3+1)/max(log(tt3+1))
}
if( tfreq == "sqrt" ){
W2trans <- tt3^(1/tfn) / max( tt3^(1/tfn) )
}
W3 <- W
X3 <- X
}else{
if( tfreq == "const" ){
W2trans <- 1*(tt3>0)
}
if( tfreq == "log" ){
W2trans <- log(tt3+1)/max(log(tt3+1))
}
if( tfreq == "sqrt" ){
W2trans <- tt3^(1/tfn) / max( tt3^(1/tfn) )
}
W3 <- spread(W2trans,ncol=ntc)
X3 <- spread(X[,seq(1,nc,ntc),drop=FALSE],ncol=ntc) + t( t(W3) * c(0:(ntc-1)) )/ncol(X)*(1-col.gap.prop)
}
# ----- ---------------- ----- #
# ----- plotting section ----- #
# ----- ---------------- ----- #
# ----- opening the basic viewport ----- #
#s0 <- 0.02 # border
#s1 <- 0.04 # yaxis
#s2 <- 0.06 # labs
#s3 <- 0.1 # expected scale
#dev.new()
#grid.rect(gp=gpar(fill="white",col="white"))
#if(exists("rmb.id",.GlobalEnv)){
#popViewport()
#grid.rect(gp=gpar(fill="white"))
if(circular) yaxis <- FALSE
#}else{
# .GlobalEnv$rmb.id <- rpois(1,lambda=10000)
#}
vp0 <- viewport(x = 0.5, y = 0.5, width = 1- 2*s0 , height = 1 - 2*s0, just = "centre", name = "vp0")
pushViewport(vp0)
vp1 <- viewport(x = 1-yaxis*s1 - s3*(!is.null(expected) ) , y = 0, width = 1- nrl*s2 - yaxis*2*s1 - s3*(!is.null(expected)), height = 1-ncl*s2, just = c("right", "bottom"), name = "vp1")
pushViewport(vp1)
grid.rect(gp=gpar(fill=bg,col=NA))
# alpha values on a log-scale
#alpha <- spread(W2,ncol=ntc)*base.alpha
#alpha <- spread(exp(2*(W2-1)),ncol=ntc)*base.alpha
if(!("alpha" %in% labels(col.opt))){
alpha <- exp(alpha.r*(W2-1))
}else{
alpha <- matrix(col.opt$alpha, ncol=ncol(W2), nrow=nrow(W2))
}
#alpha[alpha < min.alpha] <- min.alpha
if(!is.null(expected)){
alpha[alpha < 1] <- 1
}
# ----- preparing the color matrix ----- #
if(is.null(expected) | res.val.only ){
if(any(c("hsv","hcl","rgb","seq","sequential","sqn","sqt","div","diverging","diverge",
"d","s","h","t","heat","ht","ter","terrain","Wijffelaars", "w", "q17") %in% col)){
num.col <- ntc0
colv <- getcolors(num.col,col,col.opt)
}else{
if( length(col < ntc0) ){
col <- rep(col,ntc0)
}
col <- sapply(col, function(z){
if(is.integer(z)){
return(scales::alpha(z,1))
}else{
return(z)
}
})
colv <- col[1:ntc0]
num.col <- ntc0^(spine | !eqwidth | (eqwidth & circular))
}
col.Mat <- spread( t(rep(colv , nc/ntc)),nrow=nr)
if(eqwidth & !res.val.only){
dm <- dim(col.Mat)
alpha <- spread(alpha, ncol = dm[2]/ncol(alpha))#num.col #ntc0
col.Mat <- mapply( function(x,y){
z <- col2rgb(x)/255
if(nchar(x) == 9){
alpha0 <- c(substr(x,8,8),substr(x,9,9))
alpha0 <- match(alpha0,c(0:9,LETTERS[1:6]))-1
if(!any(is.na(alpha0))){
alpha0 <- (alpha0[1]*15 + alpha0[2])/240
}else{
simpleWarning("could not identify the colors alpha value")
alpha0 <- 1
}
}else{
alpha0 <- 1
}
rgb( red = z[[1]], green = z[[2]], blue = z[[3]], alpha = y*alpha0 )
}, x = col.Mat, y = alpha)
dim(col.Mat) <- dm
}
col.Mat2 <- col.Mat
}else{
resid.mat[is.nan(resid.mat)] <- 0
resid.mat[is.na(resid.mat)] <- 0
if(is.null(resid.max)){
resid.max <- max(abs(resid.mat),na.rm=T)
}
resid.mat[resid.mat < -resid.max] <- -resid.max
resid.mat[resid.mat > resid.max] <- resid.max
rmax <- ifelse(abs(resid.max) < 10^{-6},1,resid.max)
#if( !use.expected.values ){
# resid.mat <- -resid.mat
#}
alf <- 0.05
col.Mat <- apply(resid.mat,2,function(x){
if(resid.type == "pearson"){
qx <- abs(x)
qx = 2*(1 - sapply( qx, pnorm ))
#print(qx)
qx[qx < alf] <- alf
val <- log( qx )/ log( alf )
#print(val)
return( hsv( h = 2/3*(sign(x) > 0), s = val, v = 0.7 ) )
}else{
qx <- abs(x/rmax)
#rgb(sign(x) < 0,0,sign(x) > 0,alpha = exp(2*(qx-1))*base.alpha ) #qx # cut.rs/2
qx <- qx - qx %% 1/(cut.rs-1)
alf <- (exp(qx - 1)-exp(-1))/(1-exp(-1))
return( rgb(sign(x) < 0,0,sign(x) > 0,alpha = alf ) )#qx # cut.rs/2 # base.alp
}
})
colv.grey <- hsv(h=1, s = 0, v = seq(0.9,0.6,-0.3/(ntc0-1))) #spine, eqwidth, circular?
grey.Mat <- spread( t(rep(colv.grey , nc/ntc)),nrow=nr)
if( both ){
col.Mat2 <- col.Mat
col.Mat3 <- grey.Mat
col.Mat4 <- grey.Mat
col.Mat2[resid.mat >= 0] <- hsv(0,0,0,alpha=0) # red remains
col.Mat[resid.mat <= 0] <- hsv(0,0,0,alpha=0) # blue remains
col.Mat4[resid.mat < 0] <- hsv(0,0,0,alpha=0) # grey for nonreds remains
col.Mat3[resid.mat > 0] <- hsv(0,0,0,alpha=0) # grey for nonblues remains
}else{
col.Mat2 <- col.Mat
# col.Mat2[H==0] <- NA
}
}
# ----- plotting itself ----- #
# >>> divided into eqwidth <- T | eqwidth <- F >>> tfreq <- Id/sqrt/log
# >>> in spine mode, the target categories will be plotted one after another (see the for-loop)
# >>> the 3 diff. draw-function add the rectangles for the relative frequencies, the weights and the background respectively
# here were the transformations
if( circular ){
P1 = 0*X
X3 = X2 + 1/ncol(W)/2*(1-col.gap.prop)
Y3 = Y2 + 1/nrow(W)/2*(1-row.gap.prop)
#W2sqrt = tt3^(1/2) / max( tt3^(1/2) )
#W3 = spread(W2sqrt,ncol=ntc)
arat <- ncol(W3)/nrow(W3)
draw(H2/nrow(H2)*(1-row.gap.prop),H2/ncol(H2)*(1-col.gap.prop) + width.cor,X2,Y2,alpha = 1, bg = bgs, border = rect)
if(prod(dim(H2)) > 1000){
ang <- ang/2
}
if(prod(dim(H2)) > 10000){
ang <- ang/2
}
if( disp.res ){
H0 <- apply(H0,1:2, function(z) max(0.001,z))
ratMat = sqrt(H00/H0)
if( any( ratMat[H0 > 0.001] > max.rat ) ){
warning(cat("residual with ratio", max(ratMat[H0 > 0.001],na.rm=TRUE), " > ",max.rat," occurred!"))
}
ratMat <- apply(ratMat,1:2, function(z) min(max.rat, z))
sc <- dim(W3)
if(eqwidth){
R1 = matrix(1,ncol = ncol(W3)*ntc0, nrow = nrow(W3)) # 0?, ncol, ???
mv <- 1
}else{
R1 = spread(W3/2*(1-gap.prop),ncol=ntc0)# /ncol(W3)
#maximal value/ratio?
mv <- 1/2*(1-gap.prop)#/max(sc) #ncol
}
R2 <- R1 * ratMat
rscf <- max(R2[ratMat <= max.rat], na.rm = TRUE) * 2/(1-gap.prop)# * max(sc) #ncol(W3)
if(rscf < 1){
rscf <- 1
}
R2[is.na(R2)] <- 1/2*(1-gap.prop)/max(sc)#ncol
R2 <- R2/rscf
R1 <- R1/rscf
R2 <- apply(R2, 1:2, function(z) min(z, mv))
R1 <- apply(R1, 1:2, function(z) min(z, mv))
for( i in 1:length(cat.ord) ){
xind = seq( cat.ord[i], nc*ntc0, ntc0 )
P2 = P1 + H0[,xind]
#exp
#gL1 <-
mapply( function(p1,p2,rad,bg,x,y){
cc = seq( p1, p2, (p2-p1)/ceiling( (p2-p1)*ang ))
grid.polygon(x = x+c(0,rad*cos(cc*2*pi)/sc[2]),y = y+c(0, rad*sin(cc*2*pi)/sc[1]),#ncol, nrow
gp=gpar(fill=bg, lwd=1.5, col = line.col))
}, p1 = P1, p2=P2, rad = R1[,xind], bg = col.Mat[,xind], x = X3, y = Y3 )
# obs
#gL2<-
mapply( function(p1,p2,rad,bg,x,y){
cc = seq( p1, p2, (p2-p1)/ceiling( (p2-p1)*ang ))
grid.polygon(x = x+c(0,rad*cos(cc*2*pi)/sc[2]),y = y+c(0, rad*sin(cc*2*pi)/sc[1]),
gp=gpar(fill=bg, col = line.col)
)}, p1 = P1, p2=P2, rad = R2[,xind], bg = col.Mat2[,xind], x = X3, y = Y3 )
#class(gL1) <- "gList"
#class(gL2) <- "gList"
#grid.draw(gL1)
#grid.draw(gL2)
if( both ){
#exp
#gL3<-
mapply( function(p1,p2,rad,bg,x,y){
cc = seq( p1, p2, (p2-p1)/ceiling( (p2-p1)*ang ))
grid.polygon(x = x+c(0,rad*cos(cc*2*pi)/sc[2]),y = y+c(0, rad*sin(cc*2*pi)/sc[1]),
gp=gpar(fill=bg, lwd=1.5, col = line.col)
)}, p1 = P1, p2=P2, rad = R1[,xind], bg = col.Mat3[,xind], x = X3, y = Y3 )
# obs
#gL4<-
mapply( function(p1,p2,rad,bg,x,y){
cc = seq( p1, p2, (p2-p1)/ceiling( (p2-p1)*ang ))
grid.polygon(x = x+c(0,rad*cos(cc*2*pi)/sc[2]),y = y+c(0, rad*sin(cc*2*pi)/sc[1]),
gp=gpar(fill=bg,col=line.col)
)}, p1 = P1, p2=P2, rad = R2[,xind], bg = col.Mat4[,xind], x = X3, y = Y3 )
#class(gL3) <- "gList"
#class(gL4) <- "gList"
#grid.draw(gL3)
#grid.draw(gL4)
}
P1 = P2
}
}else{
#H0 <- apply(H0,1:2, function(z) max(0.0001,z))
e1<-environment()
wedge.id <- 0
wc <- list()
for( i in 1:length(cat.ord) ){
xind = seq( cat.ord[i], nc*ntc0, ntc0 )
P2 = P1 + H0[,xind]
if( eqwidth ){
R = H2*(1-gap.prop)/2 #nrow
sc <- dim(H2)
}else{
R = W3/2*(1-gap.prop) #ncol
sc <- dim(W3)
}
mapply( function(p1,p2,rad,bg,x,y){
p1 <- min(p1,1)
p2 <- min(p2,1)
p1 <- floor(10000*p1)/10000
p2 <- floor(10000*p2)/10000
if(p1 == p2){
invisible(TRUE)
}else{
cc = seq( p1, p2, (p2-p1)/ceiling( (p2-p1)*ang ))
e1$wedge.id <- e1$wedge.id+1
if( abs(p2-p1) == 1 ){
#grid.polygon(x = x+c(rad*cos(cc*2*pi)/sc[2]),y = y+c( rad*sin(cc*2*pi)/sc[1]),
#gp=gpar(fill=bg, lwd=1.2,col=line.col))
xv <- x+c(rad*cos(cc*2*pi)/sc[2])
yv <- y+c( rad*sin(cc*2*pi)/sc[1])
sp <- length(xv)
idv <- rep(e1$wedge.id,sp)
bgv <- rep(bg,sp)[1:sp]
}else{
#grid.polygon(x = x+c(0,rad*cos(cc*2*pi)/sc[2]),y = y+c(0, rad*sin(cc*2*pi)/sc[1]),
#gp=gpar(fill=bg, lwd=1.2,col=line.col))
xv <- x+c(0,rad*cos(cc*2*pi)/sc[2])
yv <- y+c(0, rad*sin(cc*2*pi)/sc[1])
sp <- length(xv)
idv <- rep(e1$wedge.id,sp)
bgv <- rep(bg,sp)[1:sp]
}
e1$wc$x <- c(e1$wc$x,xv)
e1$wc$y <- c(e1$wc$y,yv)
e1$wc$bg <- c(e1$wc$bg,bgv)
e1$wc$id <- c(e1$wc$id,idv)
return(invisible(TRUE))
}
}, p1 = P1, p2=P2, rad = R, bg = col.Mat[,xind], x = X3, y = Y3 )
P1 <- P2
}
#print(summary(as.data.frame(wc)))
#print(table(wc$id,wc$bg))
ubg <- unique(wc$bg)
for(i in ubg){
ii <- which(wc$bg == i)
grid.polygon(wc$x[ii], wc$y[ii],id=wc$id[ii],gp=gpar(fill=wc$bg[ii], lwd=1.2,col=line.col))
}
}
}else{ # not circular
if( disp.res ){
H0 <- apply(H0,1:2, function(z) max(0.001,z))
ratMat <- sqrt(H00/H0)
ratMat <- apply(ratMat,1:2, function(z) min(max.rat, z))
}
if( disp.res & spine ){
if(eqwidth){
R1 = matrix(1,ncol = ncol(W3)*ntc0/ntc, nrow = nrow(W3))/ncol(W3)/2*(1-gap.prop) # 0?
}else{
R1 = spread(W3/ncol(W3)/2*(1-gap.prop),ncol=ntc0/ntc)
}
R2 <- R1 * ratMat
rscf <- max(R2[!is.na(R2)]) * (ncol(W3)*2/(1-gap.prop))
rscf <- max(rscf,1)
}else{
rscf <- 1
}
#R2[is.na(R2)] <- 1/ncol(W3)/2*(1-gap.prop)
#R2 <- R2/rscf
#R1 <- R1/rscf
col2Mat <- matrix(col2,nrow(H2),ncol(H2))
col2Mat[H2==0] <- NA
col2Mat[W2trans==0] <- NA
draw(H2/nrow(H2)*(1-row.gap.prop),( W2trans/ncol(W2trans)*(1-col.gap.prop) + (round(W2trans,digits=7) > 0)*width.cor )/rscf,X2,Y2,alpha = 1, bg = col2Mat, border = ifelse(eqwidth, NA, line.col))
draw(H2/nrow(H2)*(1-row.gap.prop),H2/ncol(H2)*(1-col.gap.prop) + width.cor,X2,Y2,alpha = 1, bg = bgs, border = rect)
if( spine ){
for( i in 1:length(cat.ord) ){
if(i > 1){
Y <- Y + H/nrow(H)*(1-row.gap.prop)
}
xind <- seq( cat.ord[i], nc*ntc0, ntc0 )
H <- as.matrix(H0[,xind])
if( nr < 2 ){
H <- t(H)
}
CM1 <- col.Mat[,xind]
CM1[H==0] <- NA
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop)/rscf,X3,Y,alpha = 1, bg = CM1, border = NA)
if( disp.res ){
#ratMat = sqrt(H00/H0)
max.rat = max(ratMat)
W3res = W3/ncol(W3)*(1-col.gap.prop)*ratMat[,xind]
# obs
CM2 <- col.Mat2[,xind]
CM2[H==0] <- NA
draw(H/nrow(H)*(1-row.gap.prop),W3res/rscf,X3,Y,alpha = 1, bg = CM2, border = NA)
if( both ){
#exp
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop)/rscf,X3,Y,alpha = 1, bg = col.Mat3[,xind], border = NA)
# obs
draw(H/nrow(H)*(1-row.gap.prop),W3res/rscf,X3,Y,alpha = 1, bg = col.Mat4[,xind], border = NA)
}
}
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop)/rscf,X3,Y,alpha = 1, bg = NA, lwd = 1.5, border = line.col)
}
Y <- Y0
}else{
H <- apply(H,c(1,2),function(x) return(min(x/max.scale,1)))
if( disp.res ){
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop),X3,Y,alpha = 1, bg = col.Mat, border = NA)
ratMat = sqrt(H00/H0)
max.rat = max(ratMat)
W3res = W3/ncol(W3)*(1-col.gap.prop)*ratMat#[,xind]
# obs
draw(H00/nrow(H00)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop),X3,Y,alpha = 1, bg = col.Mat2, border = NA)
if( both ){
#exp
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop),X3,Y,alpha = 1, bg = col.Mat3, border = NA)
# obs
draw(H00/nrow(H00)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop),X3,Y,alpha = 1, bg = col.Mat4, border = NA)
}
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop),X3,Y,alpha = 1, bg = NA, lwd =1.5, border = line.col)
}else{
col.Mat[H==0] <- NA
draw(H/nrow(H)*(1-row.gap.prop),W3/ncol(W3)*(1-col.gap.prop),X3,Y,alpha = 1, bg = col.Mat, border = line.col)
}#border = line.col ???
}
}
if(yaxis & !circular){
if(spine){
rscf <- 1
}
sapply(y,function(z){
grid.yaxis(at =seq(z,z+1/length(y)*(1-row.gap.prop),1/length(y)*(1-row.gap.prop)/5), label = round(seq(0,max.scale,max.scale/5)/rscf,digits = 3), main = TRUE,
edits = NULL, name = NULL,
gp = gpar(fontsize = 10*lab.cex), draw = TRUE, vp = NULL)
})
sapply(y,function(z){
grid.yaxis(at =seq(z,z+1/length(y)*(1-row.gap.prop),1/length(y)*(1-row.gap.prop)/5), label = round(seq(0,max.scale,max.scale/5)/rscf,digits = 3), main = FALSE,
edits = NULL, name = NULL,
gp = gpar(fontsize=10*lab.cex), draw = TRUE, vp = NULL)
})
}
upViewport()
# ----- ---------------- ----- #
# ----- labeling section ----- #
# ----- ---------------- ----- #
if(lab){
lab.cex <- rep(lab.cex,nv)[1:nv]
if(any(col.vars*label)){
# ----- labeling the x-axis ----- #
vpX <- viewport(x = 1-yaxis*s1 - s3*(!is.null(expected)), y = 1-ncl*s2, width = 1-nrl*s2 - yaxis*2*s1 - s3*(!is.null(expected)), height = ncl*s2, just = c("right", "bottom"), name = "vpX")
pushViewport(vpX)
cur <- 1
zc <- cumprod(col.nlvl)[ which(label[which(col.vars)]) ]
for(i in 1:ncl){
vpXX <- viewport(x = 1, y = 1 - i*1/ncl, width = 1, height = 1/ncl, just = c("right", "bottom"), name = "vpXX")
pushViewport(vpXX)
if( (i > ncl-1) & num.mode.x ){
min.tick <- floor(1000*min(data[,ind[ cind[i]] ],na.rm=T ))/1000
max.tick <- floor(1000*max(data[,ind[ cind[i]] ],na.rm=T ))/1000
#vpXXb <- viewport(x = 0.5, y = 0, width = 1, height = 0.1, just = "centre", name = "vpXXb")
#pushViewport(vpXXb)
my.grid.axis(x0=0,y0=0.1,len=1,ticks= signif(seq(min.tick,max.tick,(max.tick - min.tick)/cut[cind[i]]),5), rot=0, lab.cex = lab.cex[cind[i]])
#grid.xaxis(at = seq(0,1,1/cut[cind[i]]), label = round(seq(min.tick,max.tick,(max.tick - min.tick)/cut[cind[i]]),digits = 3), main = FALSE,
#edits = NULL, name = NULL,
#gp = gpar(), draw = TRUE, vp = NULL)
#popViewport()
grid.text( names(dset)[cind[i]],x = 0.5 , y = 5/6, just = "centre",gp=gpar(fontface = "bold",cex=1.1*lab.cex))
}else{
if( varnames ){
grid.text( names(dset)[cind[i]],x = 0.5 , y = 5/6, just = "centre",gp=gpar(fontface = "bold",cex=1.1*lab.cex))
}
#z <- (nlvl[cind])[i]*cur
z <- zc[i]
grid.text( rep( clabs[[i]] ,cur), x = X[1,seq(1,nc,nc/z)] + (X[1,seq(nc/z,nc,nc/z)] - X[1,seq(1,nc,nc/z)] +1/nc*(1-col.gap.prop))/2, y = 1/3, just = "centre",gp=gpar(cex=lab.cex[ cind[i] ]))
if( boxes ){
grid.rect( y = 1/3 , x = X[1,seq(1,nc,nc/z)] , just = c("left","centre"), width = X[1,seq(nc/z,nc,nc/z)] - X[1,seq(1,nc,nc/z)] +1/nc*(1-col.gap.prop), height = 1/2, gp = gpar(fill="transparent"))
}
cur <- z
}
popViewport()
}
upViewport()
}
if(any( (!col.vars)*label)){
# ----- labeling the y-axis ----- #
vpY <- viewport(x = nrl*s2, y = 0, width = nrl*s2, height = 1-ncl*s2, just = c("right", "bottom"), name = "vpY")
pushViewport(vpY)
cur <- 1
zr <- cumprod(row.nlvl)[ which(label[which(!col.vars)]) ]
for(i in 1:nrl){
vpYY <- viewport(x = 0 + (i-1)*1/nrl, y = 0, width = 1/nrl, height = 1, just = c("left", "bottom"), name = "vpYY")
pushViewport(vpYY)
if( (i > nrl-1) & num.mode.y ){
min.tick <- floor(1000*min(data[,ind[ rind[i]] ],na.rm=T ))/1000
max.tick <- floor(1000*max(data[,ind[ rind[i]] ],na.rm=T ))/1000
#vpYYb <- viewport(x = 1, y = 0.5, width = 0.1, height = 1, just = "centre", name = "vpYYb")
#pushViewport(vpYYb)
#grid.yaxis(at = seq(0,1,1/cut[rind[i]]), label = round(seq(min.tick,max.tick,(max.tick - min.tick)/cut[rind[i]]),digits = 3), main = TRUE,
#edits = NULL, name = NULL,
#gp = gpar(), draw = TRUE, vp = NULL)
my.grid.axis(x0=0.9,y0=0,len=1,ticks=signif(seq(min.tick,max.tick,(max.tick - min.tick)/cut[rind[i]]),5), rot=90, lab.cex = lab.cex[rind[i]])
#popViewport()
grid.text( names(dset)[rind[i]],x = 1/6 , y = 0.5, just = "centre",rot=90,gp=gpar(fontface = "bold",cex=1.1*lab.cex))
}else{
if( varnames & (nr > 1) ){
grid.text( names(dset)[rind[i]],x = 1/6 , y = 0.5, just = "centre",rot=90,gp=gpar(fontface = "bold",cex=1.1*lab.cex))
}
#z <- (nlvl[rind])[i]*cur
z <- zr[i]
grid.text( rep( rlabs[[i]] ,cur), y = Y[seq(1,nr,nr/z),1] + (Y[seq(nr/z,nr,nr/z),1] - Y[seq(1,nr,nr/z),1] +1/nr*(1-row.gap.prop))/2, x = 2/3, just = "centre",rot=90,gp=gpar(cex=lab.cex[ rind[i] ]))
if( boxes ){
grid.rect( x = 2/3 , y = Y[seq(1,nr,nr/z),1] , just = c("centre","bottom"), height = Y[seq(nr/z,nr,nr/z),1] - Y[seq(1,nr,nr/z),1] +1/nr*(1-row.gap.prop), width = 1/2, gp = gpar(fill="transparent"))
}
cur <- z
}
popViewport()
}
upViewport()
}
}
# ----- scale for expected mode ----- #
if(!is.null(expected) & !res.val.only){
vpS <- viewport(x = 1, y = 0 , width = s3, height = 0.8, just = c("right","bottom"), name = "vpS")
pushViewport(vpS)
#alpha values for the residuals on a log-scale
#alfav <- exp(seq(0,-cut.rs/2+0.5,-0.5))
if( resid.type == "pearson" ){
sval <- c(log( (1 - seq(1-alf, alf, -(1-2*alf)/(cut.rs-1) ) ) )/ log( alf ))
sval <- c(sval,0, rev(sval))
alfy <- qnorm(1-alf/2)
minv <- min(-alfy,min(resid.mat[resid.mat < 0]))
maxv <- max(alfy,max(resid.mat[resid.mat > 0]))
ys <- c( minv, seq(-alfy,alfy,2*alfy/(2*cut.rs-1)), maxv)
colv <- hsv( c(rep(0,cut.rs),0,rep(2/3,cut.rs)), sval, 0.8)
grid.rect( x = 0.1 , y = ((ys - minv)/(maxv - minv))[-length(ys)], just = c("left","bottom"), height = diff(ys)/(maxv-minv), width = 0.3 , gp = gpar(fill=colv))
grid.text( label=round(ys,digits=2), y = (ys - minv)/(maxv - minv), x = 0.45, just = c("left","centre"),rot=0,gp=gpar(cex=lab.cex))
}else{
alfav <- (exp(seq(0,-1,-1/(cut.rs-1)))-exp(-1))/(1-exp(-1))
colv <- rgb(rep(c(1,0),each=cut.rs),0,rep(c(0,1),each=cut.rs),alpha = c(alfav,rev(alfav)))#*base.alpha
grid.rect( x = 0.1 , y = seq(0,1-1/2/cut.rs,1/2/cut.rs) , just = c("left","bottom"), height = 1/2/cut.rs, width = 0.3 , gp = gpar(fill=colv))
grid.text( label=round(seq(-resid.max,resid.max,resid.max/min(10,cut.rs)),digits=2), y = seq(0,1,1/2/min(10,cut.rs)), x = 0.45, just = c("left","centre"),rot=0,gp=gpar(cex=lab.cex))
}
#grid.rect( x = 0.1 , y = 0, just = c("left","bottom"), height = 1, width = 0.3 , gp = gpar(fill=hsv(0,0,0,alpha=0.1)))
upViewport()
vpS2 <- viewport(x = 1, y = 0.8 , width = 0.1, height = 0.2 , just = c("right","bottom"), name = "vpS2")
pushViewport(vpS2)
grid.text( label=resid.type, y = 0.5, x = 0.2, just = "centre",rot=90,gp=gpar(cex=lab.cex*1.2))
grid.text( label="residuals", y = 0.5, x = 0.4, just = "centre",rot=90,gp=gpar(cex=lab.cex*1.2))
grid.text( label=paste("p.value =",format(p.value,digits=2)), y = 0.5, x = 0.7, just = "centre",rot=90,gp=gpar(cex=lab.cex*1.1))
upViewport()
}
#op <- options()
#options(show.error.messages = FALSE)
#try( upViewport(), silent = TRUE)
upViewport()
if(!is.null(vp)){
try( upViewport(), silent = TRUE )
}
#options(op)
return(invisible(TRUE))
}
# ------------------------------------------------------------------------------------------------------------------ #
# ------------------------------------------------------------------------------------------------------------------ #
space <- function(nlvl,col.vars,gap.mult = 1.1,gap.prop = 0.1,last.col.zero = F, last.row.zero = F){
# >>> function to compute the gaps widths for the rmb plot
# >>> might be extended in further versions of the package
# ----- parameter check ----- #
if( !( (gap.prop > 0)&(gap.prop < 1) ) ){
stop(simpleError("Wrong gap.prop specification!"))
}
col.nlvl <- nlvl[col.vars]
row.nlvl <- nlvl[!col.vars]
nc <- prod(col.nlvl)
nr <- prod(row.nlvl)
ncv <- length(col.nlvl)
nrv <- length(row.nlvl)
col.spaces <- rep(gap.prop,nc-1)
row.spaces <- rep(gap.prop,nr-1)
if( ncv > 1 ){
cind <- 1:(nc-1)
cur <- 1
for( i in (ncv-1):1 ){
cur <- cur*col.nlvl[i+1]
curind <- which( (cind %% cur) == 0)
col.spaces[ curind ] <- col.spaces[ curind ] * gap.mult
}
}
if( nrv > 1 ){
rind <- 1:(nr-1)
cur<-1
for( i in (nrv-1):1 ){
cur <- cur*row.nlvl[i+1]
curind <- which( (rind %% cur) == 0)
row.spaces[ curind ] <- row.spaces[ curind ] * gap.mult
}
}
if(last.col.zero){
col.spaces[which(col.spaces == gap.prop)] <- 0
col.spaces <- col.spaces/gap.prop
}
if(last.row.zero){
row.spaces[which(row.spaces == gap.prop)] <- 0
row.spaces <- row.spaces/gap.prop
}
return( list(col.spaces,row.spaces))
}
# ------------------------------------------------------------------------------------------------------------------ #
# ------------------------------------------------------------------------------------------------------------------ #
spread <- function(M,ncol = 1,nrow = 1){
# >>> this function turns each matrix entry into a nrow x ncol - matrix
M <- as.matrix(M)
M2 <- apply(M,2,function(x) rep(x,each = nrow) )
if( (nrow(M) == 1)&(nrow == 1) ){
M2 <- matrix(M2,ncol = ncol(M))
#dim(M2) <- dim(M)
}
M3 <- t(apply(t(M2),2,function(x) rep(x,each = ncol) ))
if( (ncol(M) == 1)&(ncol == 1) ){
M3 <- matrix(M3,nrow = nrow(M2))
#dim(M3) <- dim(M2)
}
return(M3)
}
prettyscale = function(breaks){
#a <- min(breaks)
#b <- max(breaks)
#n <- length(breaks)
#ep <- floor(log(abs(breaks))/log(10))
ep <- sapply(breaks, function(z){
if( abs(z) < 1 ){
floor( log(abs(1/z))/log(10) )
}else{
floor( log(abs(z))/log(10) )
}
})
ep[abs(breaks)<1] <- - ep[abs(breaks)<1]
esigns = c("+","-")[(sign(ep)<0)+1]
breaks2 <- breaks
small <- which(abs(ep) < 3)
big <- which(abs(ep) > 2)
a <- format(round(breaks2[small], 2),nsmall=2)
b <- paste(round(breaks2[big]/(10^ep[big]), digits=1),paste("e",esigns[big],format(paste(0,abs(ep)[big],sep="")),sep=""),sep="")
breaks2[big] <- b
breaks2[small] <- a
return(breaks2)
}
my.grid.axis = function(x0,y0,rot, ticks, len = 1, ltm = 1/30, clockwise = FALSE, lab.cex = 0.8, keep=7, col.axis = 1, trot = 320){
grid.lines(x=c(x0, x0+len*cos(rot/180*pi)), y = c(y0, y0+len*sin(rot/180*pi)),gp=gpar(col=col.axis))
n <- length(ticks)-1
sgn <- ifelse(clockwise,-1,1)
x.ticks.1 <- x0+seq(0,1,1/n)*len*cos(rot/180*pi)
y.ticks.1 <- y0+seq(0,1,1/n)*len*sin(rot/180*pi)
x.ticks.2 <- x.ticks.1 + sgn*ltm*len*cos(rot/180*pi+pi/2)
y.ticks.2 <- y.ticks.1 + sgn*ltm*len*sin(rot/180*pi+pi/2)
ss <- seq(1,n,n/floor(keep*4 -3))
mapply( function(x0,y0,x1,y1){
grid.lines(c(x0,x1),c(y0,y1),gp=gpar(col=col.axis))
},x0 = x.ticks.1[ss], x1 = x.ticks.2[ss], y0= y.ticks.1[ss], y1 = y.ticks.2[ss])
if(rot <= 135 & rot > 45) just = ifelse(clockwise,"left","right")
if(rot <= 225 & rot > 135) just = ifelse(clockwise,"bottom","top")
if(rot <= 315 & rot > 225) just = ifelse(clockwise,"right","left")
if(rot <= 45 | rot > 315) just = ifelse(clockwise,"left","right")#top/bottom
# ii <- round( (n-1)/4 )
# if(ii > 0){
# ids <- which( floor(1:(n-2)+ii/2) %% ii != 0)
# ticks[ids+1] = ""
# }
# tmp <- n / c(5,6)
# tick.step <- trunc(tmp)[which.min( tmp - trunc(tmp) )]
if(n > keep){
keep <- round(seq(1,n,n/(keep-1)))
ticks[-c(1,n+1,keep)] = ""
}
grid.text(ticks, x.ticks.2, y.ticks.2, just = just, rot = trot, gp=gpar(fontsize=lab.cex*10,col=col.axis))
}
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.