-
Notifications
You must be signed in to change notification settings - Fork 1
/
boxcoxTransformation.R
132 lines (120 loc) · 4.14 KB
/
boxcoxTransformation.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#' Box-Cox Transformation for Simulation Summary Values
#'
#' \code{boxcoxTransformation} Box-Cox transforms summary values from a single simulation, while
#' \code{boxcoxTransformationMatrix} Box-Cox transforms summary values from multiple simulations. The
#' output of \code{boxcoxTransformationMatrix} needs to be provided as input for \code{boxcoxTransformation}.
#'
#'
#' @param summaryValuesMatrix Matrix of summary statistics from simulations
#' @param summaryValuesVector Vector of summary statistics from a single simulation
#' @param boxcoxAddition Vector of boxcox additions from \code{boxcoxTransformationMatrix}
#' @param boxcoxLambda Vector of boxcox lambda values from \code{boxcoxTransformationMatrix}
#' @return \code{boxcoxTransformation} returns a vector of Box-Cox transformed summary statistics from a
#' single observation. \code{boxcoxTransformationMatrix} returns a matrix of Box-Cox transformed summary statistics with the
#' same dimensions as \code{summaryValuesMatrix}.
#' @author Brian O'Meara and Barb Banbury
# @references O'Meara and Banbury, unpublished; Bates et al. 2009
# @keywords boxcoxEstimation Box-Cox
#' @examples
#' \donttest{
#' set.seed(1)
#' data(simRunExample)
#'
#' # example simulation
#'
#' simDataParallel <- parallelSimulateWithPriors(
#' nrepSim = 5,
#' multicore = FALSE,
#' coreLimit = 1,
#' phy = simPhyExample,
#' intrinsicFn = brownianIntrinsic,
#' extrinsicFn = nullExtrinsic,
#' startingPriorsFns = "normal",
#' startingPriorsValues = list(c(
#' mean(simCharExample[, 1]), sd(simCharExample[, 1]))),
#' intrinsicPriorsFns = c("exponential"),
#' intrinsicPriorsValues = list(10),
#' extrinsicPriorsFns = c("fixed"),
#' extrinsicPriorsValues = list(0),
#' generation.time = 100000,
#' checkpointFile = NULL,
#' checkpointFreq = 24,
#' verbose = FALSE,
#' freevector = NULL, taxonDF = NULL)
#'
#' nParFree <- sum(attr(simDataParallel, "freevector"))
#'
#' # separate the simulation results:
#' # 'true' generating parameter values from the summary values
#' summaryValuesMat <- simDataParallel[, -1:-nParFree]
#'
#' boxTranMat <- boxcoxTransformationMatrix(
#' summaryValuesMatrix = summaryValuesMat
#' )
#' boxTranMat
#'
#' boxcoxTransformation(
#' summaryValuesVector = summaryValuesMat[, 1],
#' boxcoxAddition = boxTranMat$boxcoxAddition,
#' boxcoxLambda = boxTranMat$boxcoxLambda
#' )
#' }
#' @name boxcoxTransformation
#' @rdname boxcoxTransformation
#' @export
boxcoxTransformation <- function(
summaryValuesVector,
boxcoxAddition,
boxcoxLambda
) {
###################################
#yes, a row of summary values
for (i in 1:length(summaryValuesVector)) {
addValue <- summaryValuesVector[i] + boxcoxAddition[i]
summaryValuesVector[i] <- addValue^boxcoxLambda[i]
}
return(summaryValuesVector)
}
#' @rdname boxcoxTransformation
#' @export
boxcoxTransformationMatrix <- function (summaryValuesMatrix) {
#
#library("car", quietly = T)
#
boxcoxLambda <- rep(NA, dim(summaryValuesMatrix)[2])
boxcoxAddition <- rep(NA, dim(summaryValuesMatrix)[2])
for (i in 1:dim(summaryValuesMatrix)[2]) {
boxcoxAddition[i] <- 0
lowValue <- min(
summaryValuesMatrix[, i]) - 4 * sd(summaryValuesMatrix[, i]
)
if (lowValue <= 0) {
boxcoxAddition[i] <- 4 * abs(lowValue)
}
summaryVM <- summaryValuesMatrix[, i] + boxcoxAddition[i]
boxcoxLambda[i] <- 1
if (sd(summaryValuesMatrix[, i]) > 0) {
# this alternative is thanks to:
# https://stackoverflow.com/questions/33999512/how-to-use-the-box-cox-power-transformation-in-r/34002020
# newLambda <- makeQuiet(as.numeric(try(car::powerTransform(
# summaryVM, method = "Nelder-Mead")$lambda
# )))
#
bc <- MASS::boxcox(variable ~ 1,
data = data.frame(variable = summaryVM),
lambda = seq(-20, 20, 1/10000),
plotit = FALSE)
newLambda <- bc$x[which(max(bc$y) == bc$y)[1]]
if (!is.na(newLambda)) {
boxcoxLambda[i] <- newLambda
}
}
summaryValuesMatrix[, i] <- summaryVM^boxcoxLambda[i]
}
res <- list(
boxcoxAddition = boxcoxAddition,
boxcoxLambda = boxcoxLambda,
boxcoxSummaryValuesMatrix = summaryValuesMatrix
)
return(res)
}