-
Notifications
You must be signed in to change notification settings - Fork 1
/
VariationScenario.R
executable file
·113 lines (104 loc) · 5.04 KB
/
VariationScenario.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
#' @title Scenarios for variation in CoDa regressions models
#'
#' @description
#' Scenarios of this type are illustrated in Dargel and Thomas-Agnan (2024).
#' They allow to evaluate how the response variable (Y) in a CoDa model would evolve under a hypothetical scenario for linear changes in one explanatory variable (X).
#' When the changing explanatory variable is compositional the term "linear" is understood with respect to the geometry of the simplex.
#'
#' @details
#' The linear scenario for X is computed with [seq()] in the scalar case and with [CoDa_seq()] in the compositional case.
#' The corresponding changes in Y are computed with the prediction formula, where we exploit the fact that only in one variable is changing.
#'
#'
#' @param object an object of class "lmCoDa"
#' @param Xvar a character indicating the name of the explanatory variable that changes
#' @param Xdir either character or numeric, to indicate the direction in which Xvar should change
#' - when character this should be one of the components of X, in which case the direction is the corresponding vertex of the simplex
#' - when numeric this argument is coerced to a unit vector in the simplex
#' - (when Xvar refers to a scalar variable this argument is ignored)
#' @param obs a numeric indicating the observation used for the scenario
#' @param inc_size a numeric indicating the distance between each point in the scenario of X
#' @param n_steps a numeric indicating the number of points in the scenario
#' @param add_opposite a logical, if `TRUE` the scenario also includes changes in the opposite direction
#' @param normalize_Xdir a logical, if `TRUE` the direction `Xdir` scaled to have an Aitchison norm of 1, allowing to interpret `inc_size` as the Aitchison distance
#' @return a data.frame containing the scenario of X and the corresponding predicted values of Y
#' @author Lukas Dargel
#' @references
#' - Dargel, Lukas and Christine Thomas-Agnan, “Pairwise share ratio interpretations of compositional regression models”, Computational Statistics & Data Analysis 195 (2024), p. 107945
#' @export
#' @examples
#'
#' # ---- model with scalar response ----
#' res <- lmCoDa(YIELD ~ PRECIPITATION + ilr(TEMPERATURES), data = head(rice_yields,20))
#' VariationScenario(res, Xvar = "TEMPERATURES", Xdir = "MEDIUM", n_steps = 5)
#' VariationScenario(res, Xvar = "PRECIPITATION", n_steps = 5)
#'
#'
#' # ---- model with compositional response ----
#' res <- lmCoDa(ilr(cbind(left, right, extreme_right)) ~
#' ilr(cbind(Age_1839, Age_4064)) +
#' ilr(cbind(Educ_BeforeHighschool, Educ_Highschool, Educ_Higher)) +
#' log(unemp_rate),
#' data = head(election))
#'
#' VariationScenario(res, Xvar ="cbind(Age_1839,Age_4064)",Xdir = "Age_1839", n_steps = 5)
#' VariationScenario(res, "log(unemp_rate)", n_steps = 5)
#'
VariationScenario <- function(
object,
Xvar,
Xdir,
obs = 1,
inc_size = .1,
n_steps = 100,
add_opposite = TRUE,
normalize_Xdir = TRUE) {
stopifnot(is(object, "lmCoDa"),
is.character(Xvar) && length(Xvar) == 1,
missing(Xdir) || is.character(Xdir) || is.numeric(Xdir),
is.numeric(obs) && isTRUE(obs >= 1) && obs <= nobs(object),
is.numeric(inc_size) && length(inc_size) == 1,
is.numeric(n_steps) && length(n_steps) == 1,
isTRUE(add_opposite) || isFALSE(add_opposite))
clo2 <- function(x) x/rowSums(x)
trSry <- object$trSry
# get X0
Xvar <- check_Xvar(Xvar, trSry, "NAME_COORD")
Dx <- trSry$D[[Xvar]]
scalar_x <- Dx == 0
Kx <- trSry$LR_BASE_K[[Xvar]]
X0 <- object$model[[Xvar]]
X0 <- if (scalar_x) X0[obs] else t(X0[obs,]) %*% t(Kx) # work with clr in compositional case
# get Y0
scalar_y <- trSry$D[[1]] == 0
Ky <- trSry$LR_BASE_K[[1]]
Y0 <- fitted(object)
Y0 <- if (scalar_y) Y0[obs] else Y0[obs,] %*% t(Ky) # work with clr in compositional case
inc_seq <- seq(from = -n_steps*add_opposite, to = n_steps)
if (scalar_x) {
# IDEA allow for "log-linear" sequences
Xcoef <- trSry$COEF_CLR[[Xvar]]
Xscenario <- inc_seq*inc_size
Yscenario <- if (scalar_y) Y0 else Y0[rep(1,length(Xscenario)),,drop=FALSE]
Yscenario <- Xscenario %*% Xcoef + Yscenario
Xscenario <- X0 + Xscenario
}
if (!scalar_x) {
Xdir <- check_Xdir(Xdir, colnames(X0), normalize_Xdir)
Xcoef <- trSry$COEF_CLR[[Xvar]]
Xscenario <- CoDa_path(rep(1/Dx, Dx), comp_direc = Xdir, step_size = inc_size, n_steps = n_steps, add_opposite = add_opposite)
Xscenario <- log(as.matrix(Xscenario)) %*% clrBase(Dx)
Yscenario <- if (scalar_y) Y0 else Y0[rep(1,nrow(Xscenario)),]
Yscenario <- Xscenario %*% Xcoef + Yscenario
Xscenario <- clo2(exp(X0[rep(1,nrow(Xscenario)),] + Xscenario))
}
if (!scalar_x) Xvar <- "X"
if (!scalar_y) Yvar <- "Y" else Yvar <- rownames(trSry)[1]
if (!scalar_y) Yscenario <- clo2(exp(Yscenario))
result <- data.frame(row.names = as.character(inc_seq))
result[[Yvar]] <- Yscenario
result[[Xvar]] <- Xscenario
attr(result, "inc_size") <- inc_size
attr(result, "obs") <- obs
return(result)
}