-
Notifications
You must be signed in to change notification settings - Fork 0
/
ANOPA-contrastProportions.R
215 lines (186 loc) · 8.6 KB
/
ANOPA-contrastProportions.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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
######################################################################################
#' @title contrastProportion: analysis of contrasts between proportions using Anscombe transform.
#'
#' @md
#'
#' @description The function 'contrastProportions()' performs contrasts analyses
#' on proportion data after an omnibus analysis has been obtained with 'anopa()'
#' according to the ANOPA framework. See \insertCite{lc23;textual}{ANOPA} for more.
#'
#' @param w An ANOPA object obtained from `anopa()` or `emProportions()`;
#'
#' @param contrasts A list that gives the weights for the contrasts to analyze.
#' The contrasts within the list can be given names to distinguish them.
#' The contrast weights must sum to zero and their cross-products must equal 0
#' as well.
#'
#' @return A table of significance of the different contrasts.
#'
#' @details `contrastProportions()` computes the _F_s for the contrasts,
#' testing the hypothesis that it equals zero.
#' The contrasts are each 1 degree of freedom, and the sum of the contrasts'
#' degrees of freedom totalize the effect-being-decomposed's degree of freedom.
#'
#' @references
#' \insertAllCited{}
#'
#' @examples
#' # Basic example using a one between-subject factor design with the data in compiled format.
#' # Ficticious data present success or failure of observation classified according
#' # to the state of residency (three levels); 175 participants have been observed in total.
#'
#' # The cells are unequal:
#' minimalBSExample
#'
#' # First, perform the omnibus analysis :
#' w <- anopa( {s;n} ~ state, minimalBSExample)
#' summary(w)
#'
#' # Compare the first two states jointly to the third, and
#' # compare the first to the second state:
#' cw <- contrastProportions( w, list(
#' contrast1 = c(1, 1, -2)/2,
#' contrast2 = c(1, -1, 0) )
#' )
#' summary(cw)
#'
# # Example using the Arrington et al. (2002) data, a 3 x 4 x 2 design involving
# # Location (3 levels), Trophism (4 levels) and Diel (2 levels).
# ArringtonEtAl2002
#
# # performs the omnibus analysis first (mandatory):
# w <- anopa( {s;n} ~ Location * Trophism * Diel, ArringtonEtAl2002)
# corrected(w)
#
# # execute the simple effect of Trophism for every levels of Diel and Location
# e <- emProportions(w, ~ Trophism | Diel * Location)
# summary(e)
#
# # For each of these sub-cells, contrast the four tropisms, first
# # by comparing the first two levels to the third (contrast1), second
# # by comparing the first to the second level (contrast2), and finally by
# # by comparing the first three to the last (contrast3) :
# f <- contrastProportions( e, list(
# contrast1 = c(1, 1, -2, 0)/2,
# contrast2 = c(1, -1, 0, 0),
# contrast3 = c(1, 1, 1, -3)/3
# )
# )
# summary(f)
#'
#'
######################################################################################
#'
#' @export contrastProportions
#
######################################################################################
contrastProportions <- function(
w = NULL,
contrasts = NULL
){
##############################################################################
## STEP 1: VALIDATION OF INPUT
##############################################################################
# 1.1: Is w of the right type
if (!("ANOPAobject" %in% class(w)))
stop("ANOPA::error(31): The argument w is not an ANOPA object. Exiting...")
# 1.2: Are the contrasts of the right length
if (min(unlist(lapply(contrasts, length))) != max(unlist(lapply(contrasts, length))))
stop("ANOPA::error(32): The constrasts have differing lengths. Exiting...")
relevantlevels = if (w$type == "ANOPAomnibus") {
prod( c(w$BSfactNlevels, w$WSfactNlevels) )
} else {
stop("ANOPA::oups(1): This part not yet done. Exiting...")
}
print(relevantlevels)
if (!(all(unlist(lapply(contrasts, length)) == relevantlevels )))
stop("ANOPA::error(33): The contrats lengths does not match the number of levels. Exiting...")
# 1.3a: Are the contrasts legitimate (a) sum to zero;
if (!(all(round(unlist(lapply(contrasts,sum)),8)==0)))
stop("ANOPA::error(34): Some of the constrats do not sum to 0. Exiting...")
# 1.3b: Are the contrasts legitimate (b) cross-product sum to zero
sums = c()
for (i in names(contrasts)) {for (j in setdiff(names(contrasts), i)) {
sums <-append(sums, round(sum(contrasts[[i]]*contrasts[[j]]) ),8) } }
if (!(all(sums == 0)))
stop("ANOPA::error(35): Some of the cross-products of contrasts do not totalize 0. Exiting...")
# 1.3c: Are the contrasts legitimate (c) all oppositions sum to 1
if (!(all(round(unlist(lapply(contrasts, \(x) sum(abs(x)))),8)==2)))
stop("ANOPA::error(36): Some of the contrasts' weigth does not equal 1 (for the positive weights) or -1 (for the negative weights). Exiting...")
# 1.4: is there an acceptable number of contrasts
if (length(contrasts) > relevantlevels-1)
stop("ANOPA::error(37): There are more contrasts defined than degrees of freedom. Exiting...")
##############################################################################
## STEP 2: Run the analysis
##############################################################################
# 2.1: identify the factors
ANOPAmessage("Not yet programmed...")
return(-99)
# 2.2: perform the analysis based on ???
analysis <- if(w$type == "ANOPAomnibus") {
## Contrasts on the full design
cst1way(w$compiledData, w$freqColumn, contrasts)
} else { # "ANOPAsimpleeffects"
## Contrasts on sub-data based on nesting factor(s)
cstMway(w$compiledData, w$freqColumn, w$subFactors, w$nestedFactors, contrasts)
}
##############################################################################
# STEP 3: Return the object
##############################################################################
# 3.1: preserve everything in an object of class ANOPAobject
res <- list(
type = "ANOPAcontrasts",
formula = as.formula(w$formula),
compiledData = w$compiledData,
freqColumn = w$freqColumn,
factColumns = w$factColumns,
nlevels = w$nlevels,
clevels = w$clevels,
# new information added
results = analysis
)
class(res) <- c("ANOPAobject", class(res) )
return( res )
}
##########################################
# #
# ██╗ ██╗███████╗██████╗ ███████╗ #
# ██║ ██║██╔════╝██╔══██╗██╔════╝ #
# ███████║█████╗ ██████╔╝█████╗ #
# ██╔══██║██╔══╝ ██╔══██╗██╔══╝ #
# ██║ ██║███████╗██║ ██║███████╗ #
# ╚═╝ ╚═╝╚══════╝╚═╝ ╚═╝╚══════╝ #
# #
##########################################
# Coding contrast effects #
# is postponed until interest is shown. #
##########################################
# ################################################################
# # Lets run the orthogonal contrasts #
# ################################################################
# there is only two possible cases:
# a) the contrasts are on the full data ==> cst1way
# b) the contrasts are nested with some factor(s) ==> cstMway
cst1way <- function(cData, ss, contrasts) {
# extract the proportions and the total
# get the F statistics for each contrast
# compute the correction factor cf
# compute the p values on the corrected G as usual
# This is it! let's put the results in a table
}
cstMway <- function(cData, ss, subfact, nestfact, contrasts) {
## run cst1way on every levels of the nesting factor(s)
# in case other factors are not named, collapse over these
ANOPAmessage("Not yet programmed...")
}
# ################################################################
# # Sub-functions to get observed and expected proportions #
# ################################################################
# the null hypothesis for the conditions implicated in the contrast
# ns is a vector
contrastNull <- function(contrast, ns) {
}
# the contrast hypothesis for the conditions implicated in the contrast
# ns is a vector
contrastObsd <- function(contrast, ns) {
}