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
|
#' @importFrom stats is.empty.model model.frame model.matrix model.offset na.omit na.pass reformulate relevel terms
Diff <- function(player1, player2, formula = NULL, id = "..", data = NULL,
separate.ability = NULL, refcat = NULL, contrasts = NULL,
subset = NULL) {
player.one <- player1[[id]]
player.two <- player2[[id]]
if (!is.factor(player.one) || !is.factor(player.two) ||
!identical(levels(player.one), levels(player.two)))
stop("'player1$", id, "' and 'player2$", id,
"' must be factors with the same levels")
if (!identical(attr(player.one, "contrasts"),
attr(player.two, "contrasts")))
stop("'player1$", id, "' and 'player2$", id,
"' must have the same contrasts attribute")
if(is.null(formula)) formula <- reformulate(id)
players <- levels(player.one)
nplayers <- nlevels(player.one)
ncontests <- length(player.one)
D <- matrix(nrow = ncontests, ncol = nplayers)
D <- col(D) == as.numeric(player.one)
D <- D - (col(D) == as.numeric(player.two))
colnames(D) <- paste(id, players, sep = "")
fixed <- nobars(formula)
X <- offset <- missing <- term.labels <- NULL
saturated <- FALSE
sep <- list()
empty <- is.null(fixed) || is.empty.model(mt <- terms(fixed))
if (!empty) {
factors <- attr(mt, "factors")
term.labels <- as.character(colnames(factors))
vars <- rownames(factors)
indexed <- grep("[[][^],]+[],]", vars)
if (length(indexed)) { #set NAs to zero
indices <- gsub("[^[]*[[]([^],]+)[],].*", "\\1", vars[indexed])
vars <- gsub("[[][^]]*[]]", "", vars[indexed])
## assumes no overlap, e.g. no age[..]:judge.gender[judge]
grp <- split(vars, indices)
for (ind in names(grp)) {
vars <- model.frame(terms(reformulate(grp[[ind]])),
data = data, na.action = na.pass)
lev <- levels(eval(as.name(ind), c(player1, data)))
as.sep <- rowSums(is.na(vars)) | lev %in% separate.ability
if (any(as.sep)) {
sep[[ind]] <- as.sep
vars[sep[[ind]], ] <- lapply(vars, function(x)
max(levels(x)[1], 0))
colnames(vars) <- gsub(".*[$[],? ?\"?([^]\"]*).*", "\\1",
grp[[ind]])
labels <- gsub("([^[$]*)[[$].*", "\\1", grp[[ind]])
for (lab in intersect(labels, grp[[ind]]))
data[lab] <- vars[lab]
for (lab in setdiff(labels, grp[[ind]]))
data[[lab]] <- vars[, labels == lab, drop = FALSE]
}
}
if (length(sep)) {
fixed <- reformulate(c(names(sep), attr(mt, "term.labels"),
rownames(attr(mt, "factors"))[
attr(mt, "offset")]))
mt <- terms(fixed)
}
}
idterm <- id %in% rownames(attr(mt, "factors"))
mf1 <- model.frame(mt, data = c(player1, data), na.action = na.pass)
if (nrow(mf1) != ncontests)
stop("Predictor variables are not of the correct length --",
"they probably need indexing in 'formula'.")
mf2 <- model.frame(mt, data = c(player2, data), na.action = na.pass)
if (idterm){
if (!is.null(refcat)) {
mf1[[id]] <- relevel(mf1[[id]], refcat)
mf2[[id]] <- relevel(mf2[[id]], refcat)
if (!is.null(contrasts)) contrasts[[id]] <- "contr.treatment"
} else {
## 'else' defined by contrasts arg/contrasts attr of id factor
## leave refcat NULL
if (is.null(contrasts) &
!is.null(attr(player.one, "contrasts"))){
contrasts <- list()
contrasts[[id]] <- attr(player.one, "contrasts")
}
}
}
offset <- model.offset(mf1)
if (!is.null(offset)) offset <- offset - model.offset(mf2)
if (length(sep)){ #create separate effect factor
recode <- function(x, keep){
lev <- levels(x)
ext <- make.unique(c(lev[keep], "nosep"))[sum(keep) + 1]
levels(x)[!keep] <- ext
relevel(x, ref = ext)
}
for (ind in names(grp)) {
mf1[ind] <- recode(mf1[[ind]], sep[[ind]])
mf2[ind] <- recode(mf2[[ind]], sep[[ind]])
}
}
X1 <- model.matrix(fixed, mf1, contrasts = contrasts)
X2 <- model.matrix(fixed, mf2, contrasts = contrasts)
X <- X1 - X2
## will need to check for saturation in each set of indexed var
## - however as only allowing (1|..) just consider player id for now
saturated <-
qr(na.omit(X))$rank == qr(na.omit(cbind(D, X)))$rank && !idterm
if (all(X[,1] == 0)) X <- X[, -1, drop = FALSE]
attr(X, "assign") <- attr(X1, "assign")[-1]
}
random <- findbars(formula[[2]])
if (!is.null(random)) {
if (!is.list(random)) random <- list(random)
if (length(random) > 1 ||
random[[1]] != parse(text = paste("1|", id, sep = ""))[[1]])
stop("Currently '(1 | ", id, ")' is the only random effects",
"structure allowed.")
random <- D
}
else if (!empty && (!idterm & !saturated))
warning("Ability modelled by predictors but no random effects",
call. = FALSE)
if (length(sep)) {
attr(X, "assign") <- attr(X, "assign") - 1
if (!is.null(random))
random <- D[,!sep[[id]], drop = FALSE]
}
list(X = X, random = random, offset = offset,
term.labels = term.labels, refcat = refcat, contrasts = contrasts,
saturated = saturated)
}
|
