1corrPlot <- function(design, scale="corr", recode=TRUE, cor.out=TRUE, mm.out=FALSE,
2     main.only = TRUE, three = FALSE, run.order=FALSE,
3     frml=as.formula(ifelse(three, ifelse(run.order, "~ run.no + .^3", "~ .^3"),
4                                   ifelse(run.order, "~ run.no + .^2", "~ .^2"))),
5     pal=NULL, col.grid="black", col.small="grey", lwd.grid=1.5, lwd.small=0.5,
6     lty.grid=1, lty.small=3, cex.y=1, cex.x=0.7, x.grid=NULL,
7     main=ifelse(scale=="corr","Plot of absolute correlations",ifelse(scale=="R2", "Plot of squared correlations",
8        "Plot of absolute correlations of coefficient estimates")),
9     split=0, ask=(split>0), ...){
10if (!"design" %in% class(design)) design <- data2design(design)
11if (!scale %in% c("corr", "R2", "corr.est")) stop("invalic choice for scale")
12hilf <- GWLP(design, k=5)
13if (!round(hilf[2],5)==0) stop("corrPlot is applicable for balanced designs only")
14if (!round(hilf[3],5)==0) res=2
15else {
16    if (!round(hilf[4],5)==0) res=3
17        else {
18             if (!round(hilf[5],5)==0) res=4
19             else {
20              if (!round(hilf[6],5)==0) res=5
21              else stop("The design has resolution larger than V. corrPlot is not applicable.")
22           }}}
23
24if (!is.logical(recode)) stop("recode must be logical")
25if (!is.logical(main.only)) stop("main.only must be logical")
26if (!is.logical(three)) stop("three must be logical")
27if (!is.null(x.grid) & !is.numeric(x.grid)) stop("x.grid must be a numeric vector of positions for vertical lines")
28
29if (recode) design <- change.contr(design, contrasts="contr.XuWu")
30if (res > 2 & !run.order) suppress <- TRUE else suppress <- FALSE  ## suppress main effects on horizontal axis
31if (run.order) run.no <- scale(1:nrow(design))*sqrt(nrow(design))
32
33
34design <- design[, names(factor.names(design))] ## only done here in order to have intact design structure before
35nfac <- ncol(design)
36
37mm <- model.matrix(frml, design)
38
39effid <- attr(mm, "assign")[-1]
40cpr <- abs(cor(mm[,-1]))
41if (scale=="corr.est"){
42   if (!ncol(mm)==qr(mm)$r) stop("scale=corr.est not possible because of rank deficiency")
43   cpr <- abs(cov2cor(solve(crossprod(mm))))[-1,-1]
44}
45if (mm.out || cor.out) aus <- cpr
46   if (mm.out) attr(aus, "mm") <- mm
47if (scale=="R2") cpr <- cpr^2
48diag(cpr) <- NA
49
50nc <- ncol(cpr)
51
52mecols <- which(effid %in% 1:nfac)
53if (three & !main.only){
54   int2cols <- which(effid %in% (nfac+1):((nfac*(nfac+1))/2))
55   intcols <- which(effid %in% (nfac+1):(nfac*(nfac+1)/2 + nfac*(nfac-1)*(nfac-2)/6 ))
56}
57else
58intcols <- setdiff(1:nc, mecols)   ## only 2-factor interactions present
59
60if (!three) int2cols <- intcols
61
62if (suppress & main.only) cpr <- cpr[rev(mecols), intcols]
63if (suppress & !main.only) cpr <- cpr[rev(c(mecols, int2cols)), intcols]
64if (!suppress & main.only) cpr <- cpr[rev(mecols), c(mecols, intcols)]
65if (!suppress & !main.only) cpr <- cpr[rev(c(mecols, int2cols)), c(mecols, intcols)]
66
67if (main.only) hsmall <- 0.5 + 0:length(mecols)
68else hsmall <- 0.5 + 0:(length(mecols) + length(int2cols))
69if (suppress) vsmall <- 0.5 + 0:length(intcols)
70else vsmall <- 0.5 + 0:nc
71
72if (main.only)
73   hbig <- 0.5 + c(0, which(!effid[rev(mecols[-1])]==effid[rev(mecols)[-1]]))
74else
75   hbig <- 0.5 + c(0, which(!effid[rev(c(mecols,int2cols)[-1])]==effid[rev(c(mecols,int2cols))[-1]]))
76
77#if (suppress)
78#   vbig <-
79
80if (is.null(pal)){
81if (requireNamespace("RColorBrewer"))
82 pal <- RColorBrewer::brewer.pal(9,"Blues")
83else pal <- c("white",rev(heat.colors(9)))
84}
85
86ns <- 1
87from <- 1
88to <- ncol(cpr)
89if (split>0){
90   if (ncol(cpr)/split > 8) stop("too many plots")
91   ns <- ncol(cpr)%/%split + 1
92   from <- to <- rep(0,ns)
93   for (i in 1:ns){
94       from[i] <- 1+(i-1)*split
95       to[i] <- i*split
96   }
97   to[ns] <- ncol(cpr)
98}
99
100if (!is.null(x.grid)) vbig <- x.grid
101else vbig <- 0
102
103op <- par(ask=ask)
104for (i in 1:ns){
105print(levelplot(t(round(cpr[,from[i]:to[i]],4)), col.regions=pal, cuts=length(pal)-1, aspect="fill",
106   scales=list(x=list(rot=90, tck=0, cex=cex.x), y=list(tck=0, cex=cex.y)),
107   panel = function(...){
108            panel.levelplot(...)
109            panel.abline(h = hsmall,col=col.small, lwd=lwd.small, lty=lty.small)
110            panel.abline(v = vsmall, col=col.small, lwd=lwd.small, lty=lty.small)
111            panel.abline(h = hbig,col=col.grid,lwd=lwd.grid,lty=lty.grid)
112            panel.abline(v = vbig,col=col.grid,lwd=lwd.grid,lty=lty.grid)
113        },
114   main=main,
115   xlab="",ylab=""))
116}
117par(op)
118
119if (mm.out || cor.out) invisible(aus)
120}
121