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
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
|
# $Id: expressions.R 485 2022-04-25 14:35:22Z thothorn $
### determine if an expression `x' can be interpreted as numeric
is_num <- function(x) {
if (length(x) == 1) return(is.numeric(x))
if (length(x) == 2) return(is.name(x[[1]]) && is.numeric(x[[2]]))
return(FALSE)
}
### expressions
as.char <- function(ex) {
if (length(ex) == 1) return(as.character(ex))
if (length(ex) == 3 && ex[[1]] == ":")
return(paste(as.char(ex[[2]]), ":",
as.char(ex[[3]]), sep = ""))
stop("multcomp::as.char: Failed to convert expression ", ex, " to character")
}
### extract left hand side of an expression
lhs <- function(ex) {
if (length(ex) != 1)
stop("multcomp:::lhs: expression is not of length 1", call. = FALSE )
if (length(ex[[1]]) != 3)
stop("multcomp:::lhs: expression ", sQuote(ex),
" does not contain a left and right hand side", call. = FALSE )
return(ex[[1]][[2]])
}
### extract right hand side of an expression
rhs <- function(ex) {
if (length(ex) != 1)
stop("multcomp:::rhs: expression is not of length 1", call. = FALSE )
if (length(ex[[1]][[3]]) == 2)
return(-ex[[1]][[3]][[2]])
rhs <- ex[[1]][[3]]
if (!is_num(rhs) || length(rhs) > 1)
stop("multcomp:::rhs: right hand side of expression ", sQuote(ex),
" is not a scalar numeric", call. = FALSE )
return(rhs)
}
### extract direction of the _alternative_
side <- function(ex) {
side <- as.char(ex[[1]][[1]])
if (!(side %in% c("<=", ">=", "==", "=")))
stop("multcomp:::side: does not contain ", sQuote("<=, >=, =="), call. = FALSE )
alternative <- switch(side,
"<=" = "greater",
">=" = "less",
"==" = "two.sided",
"=" = "two.sided")
return(alternative)
}
expression2coef <- function(ex, vars, debug = FALSE) {
### uses walkCode and makeCodeWalker from codetools
m.rhs <- rhs(ex)
m.lhs <- lhs(ex)
# get_coef_attr and set_coef_attr are replacements for attr and attr<-
# expression2coef originally added 'coef' attributes to symbols directly,
# but this is no longer allowed in R. Adding/accessing 'coef' attributes
# on symbols has a "global" effect within a single invocation of
# expression2coef - they are held in an environment named symcoef.
# Previously adding/accessing 'coef' attributes on symbols had a global
# effect in the whole R session. expression2coef never allows a binary
# operation that would have the same effect (variable) in both operands,
# but if that was ever supported, the scope of the 'coef' attributes on
# symbols may have to be revisited.
symcoef <- new.env(parent = emptyenv())
get_coef_attr <- function(x) {
if (is.symbol(x))
get0(as.character(x), envir = symcoef, ifnotfound = NULL)
else
attr(x, 'coef')
}
set_coef_attr <- function(x, val) {
if (is.symbol(x)) {
if (is.null(x))
rm(as.character(x), envir = symcoef)
else
assign(as.character(x), val, envir = symcoef)
} else
attr(x, 'coef') <- val
x
}
set_coef_attr( m.lhs, 1 )
if ( debug ) {
message('expression2coef',': lhs is ', sQuote(paste0(deparse(m.lhs),collapse='')))
message('expression2coef',': rhs is ', sQuote(paste0(deparse(m.rhs),collapse='')))
}
effects <-
walkCode(m.lhs,
makeCodeWalker( # dispatch operators
handler = function(v, w) {
if ( debug ) w$trace('handler',v,w)
switch( v,
'-' = w$sub,
'+' = w$add,
'*' = w$mul,
'/' = w$div,
'(' = w$exp,
':' = w$ita,
# w$fatal('handler','operator ', sQuote(v), ' is not supported')
w$eval
)
},
is.effect<-function(x) {
# vars is visible as a parameter to the enclosing function expression2coef
as.character(x) %in% vars
},
eval = function(v,w) {
if ( debug ) w$trace('eval',v,w)
parms <- c()
for ( e in as.list(v)[-1] ) {
coef <- 1
if ( debug )
message('eval',': walking ', sQuote(e), ' with coef = ', coef)
parms <- c(parms, p <- walkCode(w$setCoef(e,coef),w) )
if ( is.effect( p ) )
w$fatal('eval','within ', sQuote(deparse(v)) , ', the term ', sQuote(p),' ',
'must not denote an effect. Apart from that, ',
'the term must evaluate to a real valued constant')
}
cparms <- c()
for ( e in parms ) {
cparms <- c(cparms,parse(text=paste(e,'*',w$getCoef(e))))
}
if ( debug ) {
dumped <- lapply(cparms, function(x,w) paste(x, 'with coef =', w$getCoef(x)),w)
message('eval',': cparms = ', w$enum(dumped))
}
res <- try ( do.call(as.character(v[[1]]), as.list(cparms)), silent=T )
if ( inherits(res, 'try-error' ))
w$fatal('eval','the evaluation of the expression ', sQuote(deparse(v)),' ',
'failed with ', dQuote(attr(res,'condition')$message) )
if ( length(res) != 1 || !is.numeric(res) || !is.finite(res) )
w$fatal('eval','the expression ', sQuote(deparse(v)),' ',
'did not evaluate to a real valued constant. ',
'Result is ', sQuote(res) )
res <- w$setCoef(res*w$getCoef(v),1)
if ( debug ) {
dumped <- lapply(res, function(x,w) paste(x, 'with coef =', w$getCoef(x)),w)
message('eval',': res = ', w$enum(dumped))
}
res
},
# call -- evaluate construct directly ( this should not be reached )
call = function(v,w) {
if ( debug ) w$trace('call',v,w)
w$fatal('call',"there is probably a syntax error within subexpression", sQuote(deparse(v)))
},
# 'a - b' or '-a' -- support for subtraction of effects or constants (but not both)
sub = function(v,w) {
if ( debug ) w$trace('sub',v,w)
uminus <- length(as.list(v)) == 2
minuend <- as.list(v)[2]
subtrahend <- as.list(v)[3]
if ( debug ) {
message('sub',': minuend is ', sQuote(minuend), ', coef = ', w$getCoef(minuend) )
message('sub',': subtrahend is ', sQuote(subtrahend), ', coef = ', w$getCoef(subtrahend) )
message('sub',': uminus is ', uminus)
}
exp.coef <- ifelse( uminus, -w$getCoef(v), w$getCoef(v) )
res <- c()
for ( e in minuend ) {
if ( debug )
message('sub',': walking minuend ', sQuote(e), ', coef = ', exp.coef)
res <- c( res, walkCode(w$setCoef(e, exp.coef ),w))
}
if ( ! uminus ) {
# if uminus becomes true, subtrahend is a list of nulls.
# As a consequence, e would become null and the program would fail
for ( e in subtrahend ) {
if ( debug )
message('sub',': walking subtrahend ', sQuote(e), ', coef = ', -exp.coef)
res <- c( res, walkCode(w$setCoef(e, -exp.coef),w))
}
}
sum <- 0
symbols <- c()
# split result set into constants and symbols
for ( e in res ) {
if ( is.numeric(e) )
sum <- sum + e
else
symbols <- c(symbols, e)
}
# do not allow a reference a single effect to occur multiple times.
# to do: could be folded into a single effect by summing up coeffs
if ( length(dups <- symbols[duplicated(symbols)]) ) {
w$fatal('sub','multiple occurence of ', w$enum(dups), ' ',
'found within expression ', sQuote(deparse(v)))
}
# fold constants into single number
if ( length(symbols) == 0 ) {
return(w$setCoef(sum,1))
}
if ( sum ) {
w$fatal('sub','forming a difference between a constant and ',
'an effect as in ', sQuote(deparse(v)), ' ',
'is not supported')
}
symbols
},
# `:` `a` `b` -- support for interaction of effects as in A:B:C:D
# `:` `-a` `b` -- also support one or more signs before the first term
ita = function(v,w) {
if ( debug ) w$trace('ita',v,w)
tmp <- deparse(v)
prefix <- gsub('^([+-]*)(.*)','\\1',tmp)
name <- gsub('^([+-]*)(.*)','\\2',tmp)
sign <- 1
if ( prefix != '' ) {
for ( x in base::unlist(base::strsplit(prefix,'')) ) {
sign <- sign * switch( x,
'-' = -1,
'+' = 1,
w$fatal('ita', 'strange character ', sQuote(x),
' seen in ', sQuote(deparse(v))))
}
}
res <- w$setCoef(as.name(name), sign*w$getCoef(v) )
if ( debug ) {
dumped <- lapply(res, function(x,w) paste(x, 'with coef =', w$getCoef(x)),w)
message('ita',': res = ', w$enum(dumped))
}
res
},
# (expression)
exp = function(v,w) {
if ( debug ) w$trace('exp',v,w)
res <- c()
for ( e in as.list(v)[-1] ) {
res <- c( res, walkCode(w$setCoef(e,1),w) )
}
symbols <- c()
for ( e in res ) {
symbols <- c( symbols, w$setCoef(e, w$getCoef(e) * w$getCoef(v) ) )
}
if ( debug ) {
dumped <- lapply(symbols, function(x,w) paste(x, 'with coef =', w$getCoef(x)),w)
message('exp',': res = ', w$enum(dumped))
}
symbols
},
# '+ a b' or '+ a' -- support for the addition of constants or effects (but not both)
add = function(v,w) {
if ( debug ) w$trace('add',v,w)
res <- c()
for ( e in as.list(v)[-1] ) {
res <- c( res, walkCode(w$setCoef(e,1),w) )
}
symbols <- c()
sum <- 0
for ( e in res ) {
if ( is.numeric(e) )
sum <- sum + e
else
symbols <- c( symbols, e )
}
# fold constants into single number
if ( length(symbols) == 0 ) {
return( w$setCoef( sum * w$getCoef(v), 1 ) )
}
# do not allow to reference a single effect multiple times.
# to do: could be folded into a single effect by summing up coeffs
if ( length(dups <- symbols[duplicated(symbols)]) != 0 ) {
w$fatal('add','multiple occurence of ', w$enum(dups),' ',
'within subexpression ', sQuote(deparse(v)))
}
if ( sum ) {
w$fatal('add','adding up a constant and an effect ',
'as in ', sQuote(deparse(v)), ' is not supported')
}
# associate expression coefficient with all leafs
res <- c()
for ( e in symbols ) {
res <- c( res, w$setCoef(e, w$getCoef(e) * w$getCoef(v) ) )
}
res
},
# '* a b' -- support multiplication of constants or multiplication of an effect by a constant
mul = function(v,w) {
if ( debug ) w$trace('mul',v,w)
# collect all leafs, including constants
res <- c()
for ( e in as.list(v)[-1] ) {
res <- c(res, walkCode(w$setCoef(e,1),w) )
}
# fold literals
product <- 1
symbols <- c()
for ( r in res ) {
if ( is.numeric(r) )
product <- product * r
else
symbols <- c(symbols,r)
}
if ( product == 0 && length(symbols) ) {
w$fatal('mul','The constant part of the expression ', sQuote(deparse(v)),' ',
'evaluates to zero. This would zero out the effect(s) ', sQuote(symbols) )
}
# also take the expression coefficient into account
product <- product * w$getCoef(v)
# if only literals, return constant folding result as single number
if ( length(symbols) == 0 ) {
return(w$setCoef(product,1))
}
# prevent multiplication of fixed effects as in 'A * B' but still
# allow the multiplication effects by a constant
if ( length(symbols) > 1 && all(unlist(lapply(res,is.symbol)) ) ) {
w$fatal('mul','the multiplication of effects ', w$enum(symbols),' ',
'as in ', sQuote(deparse(v)), ' is not supported')
}
# associate the folded real valued literals as a coefficient with all symbols
res <- c()
for ( s in symbols ) {
res <- c(res, w$setCoef(s, w$getCoef(s) * product ))
}
if ( debug ) {
dumped <- lapply(res, function(x,w) paste(x, 'with coef =', w$getCoef(x)),w)
message('mul',': res = ', w$enum(dumped))
}
res
},
# '/ a b' -- division of an expression 'a' by a constant expression 'b'
div = function(v,w) {
if ( debug ) w$trace('div',v,w)
# const / const is allowed
# fixed / const is allowed: -> coef(fixed) <- 1/const
# const / fixed is forbidden
# collect all leafs, including constants
if ( (lv<-length(v)) != 3 ) {
w$fatal('div', 'internal error: length of language object ', sQuote(v), ' ',
'is not 3, but ', lv,'. Please file a bug report')
}
dividend <- c()
for ( e in as.list(v)[2] ) {
dividend <- c(dividend, walkCode(w$setCoef(e,1),w))
}
divisor <- c()
for ( e in as.list(v)[3] ) {
divisor <- c(divisor, walkCode(w$setCoef(e,1),w))
}
if ( length(divisor) != 1 ) {
w$fatal('div', "can't divide by ", sQuote(divisor), ' in ', sQuote(deparse(v)))
}
if ( any(unlist(lapply(divisor,is.effect))) ) {
w$fatal('div', "cant't divide by effect ", sQuote(divisor), ' in ', sQuote(deparse(v)))
}
if ( any(unlist(lapply(divisor,is.symbol))) ) {
w$fatal('div', "cant't divide by symbol ", sQuote(divisor), ' in ', sQuote(deparse(v)))
}
divisor <- as.numeric(divisor)
if ( !is.finite(divisor) || divisor == 0) {
w$fatal('div', "can't divide by ", sQuote(divisor), ' in ', sQuote(deparse(v)))
}
res <- c()
for ( s in dividend ) {
if ( is.numeric(s) ) {
res <- c(res, s * w$getCoef(v) / divisor )
} else {
res <- c(res, w$setCoef(s, w$getCoef(s) * w$getCoef(v) / divisor ) )
}
}
if ( debug ) {
message("div",": dividend = ", w$enum(dividend))
message("div",": divisor = ", w$enum(divisor))
message("div",": res = ", w$enum(res))
}
res
},
# leaf(e,w) -- gets called with `e` being either an effect name or a literal
leaf = function(e, w) {
if ( debug ) w$trace('leaf',e,w)
# leafs holding real valued constants tend to lose the coefficient
# attribute during implicit conversions. Hence, multiply the coefficient
# with the value and set the coefficient to one.
if ( is.numeric(e) ) {
return(w$setCoef(e * w$getCoef(e),1))
}
e
},
# return associated coefficient, or 1 if coefficient was not set before
getCoef = function(e) {
a <- get_coef_attr(e)
ifelse( is.null(a), 1, a )
},
# set coefficient of `e` to `coef` and return e
setCoef = function(e,coef) {
set_coef_attr(e, coef)
},
enum = function(x) {
paste0("'" ,x, "'",collapse=', ')
},
fatal = function(name,...) {
stop(paste0('multcomp:::expression2coef::walkCode::',name),': ', ..., call. = FALSE )
},
trace = function(fn,v,w) {
message(fn,': v = ', sQuote(v),
', mode = ', mode(v),
', typeof = ', typeof(v),
', length = ', length(v),
', coef = ', w$getCoef(v))
}
)) # end of walkCode(lhs, makeCodeWalker( ... ) )
if ( any(idx <- is.numeric(effects) ) ) {
stop('multcomp:::expression2coef: The lhs expression ', sQuote(deparse(m.lhs)), ' ',
'contains a numeric offset term evaluating to ', paste0(effects[idx],collapse=', '), '. ',
'This is either an internal error or a misspecification from your part. ',
'If so, please pull these offsets to the right-hand side of the equation', call. = FALSE )
}
effect.names <- c()
effect.coefs <- c()
# There might be only a single effect as in 'Agriculture = 0'. Thus use
# c(effects) to prevent the for loop from running into an error condition
for ( effect in c(effects) ) {
effect.names <- c( effect.names, as.character(effect))
effect.coefs <- c( effect.coefs, get_coef_attr(effect))
}
list( coef = effect.coefs,
names = effect.names,
m = m.rhs,
alternative = side(ex),
lhs = deparse( m.lhs, width.cutoff = 500 ) )
}
### interpret character representations of linear functions
chrlinfct2matrix <- function(ex, var) {
if (!is.character(ex))
stop("multcomp:::chrlinfct2matrix: argument ", sQuote(ex), " is not of type character", call. = FALSE )
if (!is.character(var))
stop("multcomp:::chrlinfct2matrix: argument ", sQuote(var), " is not of type character", call. = FALSE )
K <- matrix(0, nrow = length(ex), ncol = length(var))
colnames(K) <- var
rownames(K) <- seq_along(ex)
m <- rep(0, length(ex))
for (i in 1:length(ex)) {
expr <- parse(text = ex[i])
if (length(expr[[1]]) != 3)
stop("multcomp:::chrlinfct2matrix: argument ", sQuote(ex[i]),
" cannot be interpreted as expression", call. = FALSE )
tmp <- expression2coef(expr,vars=var)
### (Intercept) lost () in expression2coef
if ("(Intercept)" %in% var)
tmp$names[tmp$names == "Intercept"] <- "(Intercept)"
if (!all(tmp$names %in% var))
stop("multcomp:::chrlinfct2matrix: variable(s) ",
paste(sQuote(tmp$names[!tmp$names %in% var]),collapse=', '), " not found", call. = FALSE )
for (n in tmp$names)
K[i, var == n] <- tmp$coef[tmp$names == n]
m[i] <- tmp$m
if (i == 1) {
alternative <- tmp$alternative
} else {
if (tmp$alternative != alternative)
stop("multcomp:::chrlinfct2matrix: mix of alternatives currently not implemented", call. = FALSE )
}
rownames(K)[i] <- paste0(tmp$lhs, collapse = "")
}
list(K = K, m = m, alternative = alternative)
}
|
