diff --git a/R/profileInformationIndex.R b/R/profileInformationIndex.R index b81765ac9..247a0edc8 100644 --- a/R/profileInformationIndex.R +++ b/R/profileInformationIndex.R @@ -1,5 +1,5 @@ -# x: vector of any type +# x: vector of any type, NA removed # d: number of significant digits to retain for numeric x .prepareVector <- function(x, d) { @@ -31,7 +31,8 @@ return(x) } - +# x: vector of any type, ignoring NA +# d: number of significant digits to retain for numeric x .prepareBaseline <- function(x, d, type) { if(type == 'numeric') { @@ -57,7 +58,7 @@ # integer-coded (factor, logical) # frequency weighted by hz thickness, # non-NA value - b <- names(sort(table(x), decreasing = TRUE))[1] + b <- names(sort(table(x, useNA = 'no'), decreasing = TRUE))[1] b <- rep(b, times = length(x)) } @@ -87,20 +88,25 @@ x <- .prepareVector(x, d = d) } + # otherwise using a single string, with NA removed + # length of compressed vector is the metric res <- length(memCompress(x, type = m)) + + # length of empty string ('') is 8 bits + res <- res - 8 + return(res) } -## TODO: think more about NA handling # i: vector, any type # numericDigits: requested precision for numeric data -# removeNA: remove NA before preparing -.prepareVariable <- function(i, numericDigits, removeNA) { +.prepareVariable <- function(i, numericDigits) { - ## TODO: what is appropriate return value when all NA? + # remove all NA + i <- as.vector(na.omit(i)) # short-circuit for all NA if(all(is.na(i))) { @@ -109,10 +115,6 @@ # baseline is mean(i) if(is.numeric(i)) { - # optionally remove NA from length(gz()) - if(removeNA) { - i <- na.omit(i) - } # baseline for numeric data: 1cm slice wt. mean b <- .prepareBaseline(i, d = numericDigits, type = 'numeric') @@ -124,11 +126,6 @@ # categorical or logical # baseline is most-frequent, non-NA value - # optionally remove NA from length(gz()) - if(removeNA) { - i <- na.omit(i) - } - # treating all categorical variables as nominal for now # baseline is the most frequent, weighted by hz thickness, non-NA value @@ -143,7 +140,7 @@ } -.PII_by_profile <- function(x, vars, baseline, method, numericDigits, removeNA, compression) { +.PII_by_profile <- function(x, vars, baseline, method, numericDigits, compression) { # diced SPC -> DF of horizon data # just variables of interest @@ -153,8 +150,7 @@ h <- lapply( h[, vars, drop = FALSE], FUN = .prepareVariable, - numericDigits = numericDigits, - removeNA = removeNA + numericDigits = numericDigits ) # extract variables @@ -219,9 +215,9 @@ #' #' @param numericDigits integer, number of significant digits to retain in numeric -> character conversion #' -#' @param removeNA remove NA before compression +#' @param padNA logical, pad depths to `max(x)`, supplied to `dice(fill = padNA)` #' -#' @param padNA pad all profiles with NA to deepest lower depth in `x` (`max(x)`) +#' @param scaleNumeric logical, `scale()` each numeric variable, causing "profile information" to vary based on other profiles in the collection #' #' @param compression character, compression method as used by [memCompress()]: 'gzip', 'bzip2', 'xz', 'none' #' @@ -293,29 +289,34 @@ #' text(x = 1:5, y = 105, labels = pi, cex = 0.85) #' mtext('Profile Information Index (bytes)', side = 1, line = -1) #' -profileInformationIndex <- function(x, vars, method = c('joint', 'individual'), baseline = FALSE, numericDigits = 8, removeNA = !padNA, padNA = FALSE, compression = 'gzip') { +profileInformationIndex <- function(x, vars, method = c('joint', 'individual'), baseline = FALSE, numericDigits = 8, padNA = FALSE, scaleNumeric = FALSE, compression = 'gzip') { # sanity check method <- match.arg(method) - # scale numeric variables to mean of 0 and SD of 1 - for(i in vars) { - if(inherits(x[[i]], 'numeric')) { - - # scale() will return NaN when SD() is NA or 0 - .sd_i <- sd(x[[i]], na.rm = TRUE) - - # threshold for very small numbers - if(!is.na(.sd_i)) { - if(.sd_i > 0.0001) { - x[[i]] <- scale(x[[i]]) + ## TODO: think more about this + # scaling means that PII depends on collection and is not absolute + if(scaleNumeric) { + # scale numeric variables to mean of 0 and SD of 1 + for(i in vars) { + if(inherits(x[[i]], 'numeric')) { + + # scale() will return NaN when SD() is NA or 0 + .sd_i <- sd(x[[i]], na.rm = TRUE) + + # threshold for very small numbers + if(!is.na(.sd_i)) { + if(.sd_i > 0.0001) { + x[[i]] <- scale(x[[i]]) + } } + + # otherwise, use as-is } - - # otherwise, use as-is } } + ## TODO: this will error / drop profiles in the presence of bad horizonation # dice() to 1cm intervals for common baseline @@ -339,7 +340,6 @@ profileInformationIndex <- function(x, vars, method = c('joint', 'individual'), baseline = baseline, method = method, numericDigits = numericDigits, - removeNA = removeNA, compression = compression ) diff --git a/misc/sandbox/profile-information-index.R b/misc/sandbox/profile-information-index.R index eeb3897bf..f413e6b4d 100644 --- a/misc/sandbox/profile-information-index.R +++ b/misc/sandbox/profile-information-index.R @@ -27,8 +27,8 @@ length(memCompress(a, type = 'gzip')) / length(memCompress(b, type = 'gzip')) aqp:::.prepareVector(1:10, d = 4) -aqp:::.prepareVariable(1:10, numericDigits = 4, removeNA = FALSE) -aqp:::.prepareVariable(c('A', 'A', 'A', 'B', 'C'), numericDigits = 4, removeNA = FALSE) +aqp:::.prepareVariable(1:10, numericDigits = 4) +aqp:::.prepareVariable(c('A', 'A', 'A', 'B', 'C'), numericDigits = 4) # note usage aqp:::.compressedLength(1:10) @@ -44,6 +44,11 @@ aqp:::.compressedLength('A') aqp:::.compressedLength(factor('A')) +# NA should not "add" information +aqp:::.compressedLength(c('A', 'B', NA)) +aqp:::.compressedLength(c('A', 'B', 'C')) +aqp:::.compressedLength(NA) + aqp:::.compressedLength(0) aqp:::.compressedLength(1000) @@ -137,7 +142,7 @@ p2 <- data.frame( p3 <- data.frame( id = 3, top = c(0, 10, 20, 30, 40, 50), bottom = c(10, 20, 30, 40, 50, 100), - p = c(1, 5, 10, 35, 6, 2), + p = c(5, 8, 10, 35, 6, 2), name = c('A1', 'A2', 'Bw', 'Bt1', 'Bt2', 'C') ) @@ -145,7 +150,7 @@ p3 <- data.frame( p4 <- data.frame( id = 4, top = c(0, 10, 20, 30, 40, 50), bottom = c(10, 20, 30, 40, 50, 100), - p = c(1, NA, NA, NA, NA, NA), + p = c(5, NA, NA, NA, NA, NA), name = c('A1', 'A2', 'Bw', 'Bt1', 'Bt2', 'C') ) @@ -194,7 +199,6 @@ mtext('Profile Complexity Ratio', side = 1, line = -0.5) - # effect of aggregation function profileInformationIndex(z, vars = vars, method = 'i', baseline = TRUE) profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE) @@ -203,14 +207,11 @@ profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE) profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE) profileInformationIndex(z, vars = vars, method = 'j', baseline = FALSE) -# effect of removing NA -profileInformationIndex(z, vars = vars, method = 'j', baseline = FALSE, removeNA = TRUE, padNA = FALSE) -profileInformationIndex(z, vars = vars, method = 'j', baseline = FALSE, removeNA = FALSE, padNA = FALSE) - - -profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE, padNA = TRUE, removeNA = FALSE) -profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE, padNA = FALSE, removeNA = FALSE) +# effect of padding depths, only when vars includes top/bottom +profileInformationIndex(z, vars = c(vars, 'top'), method = 'j', baseline = FALSE, padNA = TRUE) +profileInformationIndex(z, vars = c(vars, 'top'), method = 'j', baseline = FALSE, padNA = FALSE) +# effect of number digits in character representation profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE, numericDigits = 1) profileInformationIndex(z, vars = vars, method = 'j', baseline = TRUE, numericDigits = 10) @@ -265,31 +266,23 @@ x <- fetchOSD(s) # vars <- c('hue', 'value', 'chroma', 'texture_class', 'cf_class', 'pH', 'pH_class', 'distinctness', 'topography') + par(mar = c(3, 0, 1, 2), mfrow = c(2, 1)) vars <- c('hue', 'value', 'chroma', 'hzname') -x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = TRUE, removeNA = FALSE) +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = TRUE, removeNA = FALSE, padNA = TRUE') +title('baseline = FALSE') -x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = FALSE, removeNA = TRUE) +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = TRUE) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = TRUE, removeNA = TRUE, padNA = FALSE') - - - -x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = FALSE, removeNA = TRUE) - -plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) -axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = FALSE, removeNA = FALSE, padNA = FALSE') - +title('baseline = TRUE') @@ -300,69 +293,65 @@ par(mar = c(3, 0, 1, 2), mfrow = c(2, 1)) vars <- c('hue', 'value', 'chroma', 'hzname') -x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = FALSE) +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('method = j, baseline = FALSE, padNA = FALSE', cex.main = 1) +title('method = j, baseline = FALSE', cex.main = 1) -x$pi <- profileInformationIndex(x, vars = vars, method = 'i', baseline = FALSE, padNA = FALSE) +x$pi <- profileInformationIndex(x, vars = vars, method = 'i', baseline = FALSE) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('method = i, baseline = FALSE, padNA = FALSE', cex.main = 1) +title('method = i, baseline = FALSE', cex.main = 1) - -x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = TRUE, removeNA = FALSE) +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = TRUE) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('method = j, baseline = FALSE, padNA = TRUE', cex.main = 1) +title('method = j, baseline = TRUE', cex.main = 1) -x$pi <- profileInformationIndex(x, vars = vars, method = 'i', baseline = FALSE, padNA = TRUE, removeNA = FALSE) +x$pi <- profileInformationIndex(x, vars = vars, method = 'i', baseline = TRUE) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('method = i, baseline = FALSE, padNA = TRUE', cex.main = 1) - -dev.off() - +title('method = i, baseline = TRUE', cex.main = 1) +# add top depth to test padNA = TRUE +vars <- c('hue', 'value', 'chroma', 'hzname', 'top') +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = TRUE) +plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) +axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) +title('method = j, baseline = FALSE, padNA = TRUE', cex.main = 1) -x$pi <- profileInformationIndex(x, vars = vars, baseline = TRUE, method = 'median') +x$pi <- profileInformationIndex(x, vars = vars, method = 'i', baseline = FALSE, padNA = TRUE) -par(mar = c(3, 0, 1, 2)) -plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE) +plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = TRUE, method = median') +title('method = i, baseline = FALSE, padNA = TRUE', cex.main = 1) -x$pi <- profileInformationIndex(x, vars = vars, baseline = FALSE, method = 'sum') +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = TRUE, padNA = TRUE) -par(mar = c(3, 0, 1, 2)) -plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE) +plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = FALSE, method = sum') - +title('method = j, baseline = TRUE, padNA = TRUE', cex.main = 1) -x$pi <- profileInformationIndex(x, vars = vars, baseline = TRUE, method = 'sum') +x$pi <- profileInformationIndex(x, vars = vars, method = 'i', baseline = TRUE, padNA = TRUE) -par(mar = c(3, 0, 1, 2)) -plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE) +plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE, max.depth = 200) axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = TRUE, method = sum') +title('method = i, baseline = TRUE, padNA = TRUE', cex.main = 1) + + +dev.off() -x$pi <- profileInformationIndex(x, vars = 'hue', baseline = TRUE, method = 'sum') -par(mar = c(3, 0, 1, 2)) -plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE) -axis(side = 1, at = 1:length(x), labels = format(x$pi, digits = 3)[order(x$pi)], cex.axis = 0.75, las = 1) -title('baseline = TRUE, method = sum') @@ -375,7 +364,7 @@ x$A <- .lab[, 2] x$B <- .lab[, 3] -x$pi <- profileInformationIndex(x, vars = c('L', 'A', 'B'), baseline = FALSE, method = 'j', padNA = FALSE, removeNA = TRUE) +x$pi <- profileInformationIndex(x, vars = c('L', 'A', 'B'), baseline = FALSE, method = 'j') par(mar = c(3, 0, 1, 2), mfrow = c(1,1)) plotSPC(x, width = 0.3, name.style = 'center-center', plot.order = order(x$pi), cex.names = 0.66, shrink = TRUE) @@ -420,10 +409,10 @@ vars <- c('hzname', 'hue', 'value', 'chroma', 'texture_class') z <- data.frame( - baseline.joint = profileInformationIndex(x, vars = vars, baseline = TRUE, method = 'j', padNA = FALSE, removeNA = TRUE), - baseline.individual = profileInformationIndex(x, vars = vars, baseline = TRUE, method = 'i', padNA = FALSE, removeNA = TRUE), - joint = profileInformationIndex(x, vars = vars, baseline = FALSE, method = 'j', padNA = FALSE, removeNA = TRUE), - individual = profileInformationIndex(x, vars = vars, baseline = FALSE, method = 'i', padNA = FALSE, removeNA = TRUE) + baseline.joint = profileInformationIndex(x, vars = vars, baseline = TRUE, method = 'j', padNA = FALSE), + baseline.individual = profileInformationIndex(x, vars = vars, baseline = TRUE, method = 'i', padNA = FALSE), + joint = profileInformationIndex(x, vars = vars, baseline = FALSE, method = 'j', padNA = FALSE), + individual = profileInformationIndex(x, vars = vars, baseline = FALSE, method = 'i', padNA = FALSE) ) cor(z) @@ -434,11 +423,21 @@ cor(z) hexplom(z, par.settings = tactile.theme(axis.text = list(cex = 0.66)), trans = log, inv = exp, xbins = 30, colramp = .cp, colorkey = FALSE, varname.cex = 0.75, varname.font = 2, main = 'Profile Information Index', xlab = '') +## interesting... very little difference between joint vs. individual + + plot(joint ~ individual, data = z, las = 1) plot(baseline.joint ~ baseline.individual, data = z, las = 1) hexbinplot(baseline.joint ~ joint, data = z, par.settings = tactile.theme(), trans = log, inv = exp, xbins = 30, colramp = .cp, colorkey = FALSE, varname.cex = 0.75, varname.font = 2, main = 'Profile Information Index') -x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = TRUE, padNA = FALSE, removeNA = TRUE) + +## Profile Complexity Ratio (baseline = TRUE): less influenced by number of horizons +## --> cor ~ 0.44 + +## Profile Complexity Index (baseline = FALSE): more influenced by number of horizons +## --> cor ~ 0.71 + +x$pi <- profileInformationIndex(x, vars = vars, method = 'j', baseline = FALSE, padNA = FALSE) x$nhz <- profileApply(x, FUN = nrow, simplify = TRUE) x$greatgroup <- factor(x$greatgroup, levels = c('palexeralfs', 'haploxeralfs', 'haploxerepts', 'xerorthents', 'haploxererts', 'endoaquolls')) @@ -504,7 +503,9 @@ site(z2)$g <- 'RW' z <- combine(z1, z2) -z$pi <- profileInformationIndex(z, vars = c('p1', 'p2', 'p3'), method = 'j', baseline = FALSE, padNA = FALSE, removeNA = TRUE) +## important: scales may be very different, scale() -> compress + +z$pi <- profileInformationIndex(z, vars = c('p1', 'p2', 'p3'), method = 'j', scaleNumeric = TRUE) z$nhz <- profileApply(z, FUN = nrow, simplify = TRUE) z$g <- factor(z$g) @@ -515,6 +516,13 @@ z$nhz par(mar = c(0, 0, 3, 2)) groupedProfilePlot(z, groups = 'g', color = 'p1') +lsp <- get('last_spc_plot', envir = aqp.env) + +o <- lsp$plot.order +.b <- z[, , .LAST, .BOTTOM] + +text(x = o, y = .b, labels = z$pi, cex = 0.66, pos = 1) + bwplot(g ~ pi, data = site(z)) @@ -540,54 +548,16 @@ z1 <- lapply( z1 <- combine(z1) # z1 <- trunc(z1, 0, min(z1)) -z1$pi <- profileInformationIndex(z1, vars = c('p1', 'p2', 'p3'), method = 'j', baseline = TRUE, padNA = FALSE, removeNA = TRUE) +z1$pi <- profileInformationIndex(z1, vars = c('p1', 'p2', 'p3'), method = 'j', scale = TRUE, baseline = TRUE, padNA = FALSE) par(mar = c(3, 0, 0, 2)) plotSPC(z1, color = 'p1', plot.order = order(z1$pi), print.id = FALSE, width = 0.35, divide.hz = FALSE) axis(side = 1, at = 1:length(z1), labels = format(z1$pi[order(z1$pi)], digits = 3), cex.axis = 0.66) -# simulate three profiles of increasing complexity -p1 <- data.frame(id = 1, top = 0, bottom = 100, p = 5) -p2 <- data.frame(id = 2, top = c(0, 10, 20, 30, 40, 50), bottom = c(10, 20, 30, 40, 50, 100), p = rep(5, times = 6)) -p3 <- data.frame(id = 3, top = c(0, 10, 20, 30, 40, 50), bottom = c(10, 20, 30, 40, 50, 100), p = c(1, 5, 10, 3, 6, 2)) - -# combine and upgrade to SPC -z <- rbind(p1, p2, p3) -depths(z) <- id ~ top + bottom - -# visual check -plotSPC(z, color = 'p') - -# compute information index several ways -profileInformationIndex(z, vars = c('p'), method = 'sum') -profileInformationIndex(z, vars = c('p'), method = 'mean') - -profileInformationIndex(z, vars = c('p'), method = 'mean', baseline = FALSE) -profileInformationIndex(z, vars = c('p'), method = 'sum', baseline = FALSE) - - - -## ... need to resolve this -# effect of profile depth -p1 <- data.frame(id = 1, top = 0, bottom = 75, p = 5) -p2 <- data.frame(id = 2, top = c(0, 10, 20, 30, 40, 50), bottom = c(10, 20, 30, 40, 50, 100), p = rep(5, times = 6)) -p3 <- data.frame(id = 3, top = c(0, 10, 20, 30, 40, 50), bottom = c(10, 20, 30, 40, 50, 150), p = c(1, 5, 10, 3, 6, 2)) - -# combine and upgrade to SPC -z <- rbind(p1, p2, p3) -depths(z) <- id ~ top + bottom - -# visual check -plotSPC(z, color = 'p') - -# compute information index several ways -profileInformationIndex(z, vars = c('p'), method = 'sum') -profileInformationIndex(z, vars = c('p'), method = 'mean') - -profileInformationIndex(z, vars = c('p'), method = 'mean', baseline = FALSE) -profileInformationIndex(z, vars = c('p'), method = 'sum', baseline = FALSE) - - +z1$pi <- profileInformationIndex(z1, vars = c('p1', 'p2', 'p3'), method = 'j', scale = TRUE, baseline = FALSE, padNA = FALSE) +par(mar = c(3, 0, 0, 2)) +plotSPC(z1, color = 'p1', plot.order = order(z1$pi), print.id = FALSE, width = 0.35, divide.hz = FALSE) +axis(side = 1, at = 1:length(z1), labels = format(z1$pi[order(z1$pi)], digits = 3), cex.axis = 0.66)