-
Notifications
You must be signed in to change notification settings - Fork 1
/
class-lmCoDa.R
executable file
·177 lines (152 loc) · 5.27 KB
/
class-lmCoDa.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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# ---- constructors -----------------------------------------------------------
#' Estimating CoDa regression models
#'
#' This is a thin wrapper around [lm()] followed by [ToSimplex()], which
#' allows to create a lmCoDa object in one step.
#'
#' @param formula as in [lm()]
#' @param data as in [lm()]
#' @param ... arguments passed on to [lm()]
#' @return an object of class "lm" and "lmCoDa" if the formula include at least
#' one log-transformation
#'
#' @author Lukas Dargel
#' @seealso [lm()], [ToSimplex()], [compositions::ilr()], [compositions::alr()]
#' @export
#' @examples
#'
#' # XY-compositional model
#' res <- lmCoDa(
#' ilr(cbind(left, right, extreme_right)) ~
#' ilr(cbind(Educ_BeforeHighschool, Educ_Highschool, Educ_Higher)),
#' data = head(election, 20))
#'
#' # X-compositional model
#' res <- lmCoDa(YIELD ~ PRECIPITATION + ilr(TEMPERATURES), data = head(rice_yields, 20))
#'
lmCoDa <- function(formula, data, ...) {
cl <- match.call()
cl[[1]] <- as.symbol("lm")
res <- eval(cl)
cl[[1]] <- as.symbol("lmCoDa")
res$call <- cl
return(ToSimplex(res))
}
#' Converting Linear Models to CoDa models
#'
#' @description
#' The function converts the output of a "lm" to the "lmCoDa" class, which
#' offers additional tools for the interpretation of a CoDa regression models.
#' Most of the work is done by the [transformationSummary()] function, which
#' has its own documentation page, but should be reserved for internal use.
#'
#'
#' @inherit lmCoDa return examples
#' @inheritParams VariationScenario
#' @seealso [lm()], [lmCoDa()]
#' @author
#' - Lukas Dargel
#' - Rodrigue Nasr
#' @export
#' @examples
#'
#' # XY-compositional model
#' res <- lm(
#' ilr(cbind(left, right, extreme_right)) ~
#' ilr(cbind(Educ_BeforeHighschool, Educ_Highschool, Educ_Higher)),
#' data = head(election, 20))
#' res <- ToSimplex(res)
#'
#' # X-compositional model
#' res <- lm(YIELD ~ PRECIPITATION + ilr(TEMPERATURES), data = head(rice_yields, 20))
#' res <- ToSimplex(res)
ToSimplex <- function(object){
trSry <- transformationSummary(object)
if (all(trSry[,'LR_TRAN'] == ""))
return(object)
object$trSry <- trSry
class(object) <- c('lmCoDa', class(object))
if (trSry[1,'LR_TRAN'] == "")
return(object)
invTran <- match.fun(paste0(trSry[1,'LR_TRAN'],"Inv"))
attributes(object$fitted.values)$orig <- invTran(object$fitted.values)
attributes(object$residuals)$orig <- invTran(object$residuals)
return(object)
}
# ---- methods ----------------------------------------------------------------
#' @inherit predict.lmCoDa title description details params
#' @param split logical, if `TRUE` the coefficients are reported as a list instead
#' of a matrix, where list structure reflects the explanatory variables of the model
#' @param ... not used
#' @return a matrix
#'
#' @author Lukas Dargel
#' @exportS3Method
coef.lmCoDa <- function(object, space = NULL, split = FALSE, ...) {
stopifnot(is.null(space) || space %in% c("clr", "simplex"))
type <- if (is.null(space)) "COEF_COORD" else c("clr" = "COEF_CLR", "simplex" = "COEF_SIMPLEX")[space]
cfs <- object$trSry[[type]][-1]
cfs <- if (split) cfs else Reduce("rbind", cfs)
return(cfs)
}
#' @inherit predict.lmCoDa title description details params
#' @return matrix or vector
#' @importFrom stats fitted
#' @author Lukas Dargel
#' @exportS3Method
fitted.lmCoDa <- function(object, space = NULL, ...) {
stopifnot(is.null(space) || space %in% c("clr", "simplex"))
fity <- object$fitted.values
if (is.null(space))
return(fity)
Ky <- object$trSry[["LR_BASE_K"]][[1]]
if (length(Ky) == 0)
stop("The space argument can only be used for Y compositional models!")
if (space == "simplex")
return(attr(fity, "orig"))
if (space == "clr")
return(fity %*% t(Ky))
}
#' Predictions, fitted values, residuals, and coefficients in CoDa models
#'
#' These functions work as in the usual lm object.
#' They additionally offer the possibility use the `space` argument
#' which transforms them into directly into clr space or in the simplex.
#'
#' @param object class "lmCoDa"
#' @param space a character indicating in which space the prediction should
#' be returned. Supported are the options `c("clr", "simplex")`.
#' @param ... passed on to [predict.lm()]
#' @return matrix or vector
#' @author Lukas Dargel
#' @exportS3Method
predict.lmCoDa <- function(object, space = NULL, ...) {
stopifnot(is.null(space) || space %in% c("clr", "simplex"))
pred <- NextMethod("predict")
if (is.null(space))
return(pred)
Ky <- object$trSry[["LR_BASE_K"]][[1]]
if (length(Ky) == 0)
stop("The space argument can only be used for Y compositional models!")
if (space == "simplex")
return(clrInv(pred %*% t(Ky)))
if (space == "clr")
return(pred %*% t(Ky))
}
#' @inherit predict.lmCoDa title description details params
#' @return matrix or vector
#' @author Lukas Dargel
#' @exportS3Method
residuals.lmCoDa <- function(object, space = NULL, ...) {
stopifnot(is.null(space) || space %in% c("clr", "simplex"))
resi <- object$residuals
if (is.null(space))
return(resi)
Ky <- object$trSry[["LR_BASE_K"]][[1]]
if (length(Ky) == 0)
stop("The space argument can only be used for Y compositional models!")
if (space == "simplex")
return(attr(resi, "orig"))
if (space == "clr")
return(resi %*% t(Ky))
}