Description Usage Arguments Examples
mudball performs automatic backward elimination of all effects of linear mixed effect model and returns a list of the backward eliminated models.
1 |
model |
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 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (model, i = 0, beeping = F, ps = list())
{
require(magrittr)
require(lme4)
require(beepr)
if (beeping) {
beep()
}
message("###########################################")
print(i)
model_summary = summary(model)
rand_factors = (model_summary$ngrps %>% as.data.frame() %>%
rownames() %>% as.list())
variances = model_summary$varcor
stds = list()
for (rand_factor in rand_factors) {
rand_factor_std = sprintf("std_%s", rand_factor)
get_rand_factor_std = sprintf("attr(variances$%s,\"stddev\")",
rand_factor)
assign(rand_factor_std, as.data.frame(eval(parse(text = get_rand_factor_std))))
stds = append(stds, rand_factor_std)
}
min_std = 10000
min_list = list()
for (rand_factor in rand_factors) {
if (nrow(get(sprintf("std_%s", rand_factor))) == 1)
next
tmp_min = min(get(sprintf("std_%s", rand_factor))[-1,
])
if (tmp_min <= min_std) {
min_std = tmp_min
list_of_vars = get(sprintf("std_%s", rand_factor))[,
1]
min_index = which(list_of_vars == min_std)[1]
min_name = rownames(get(sprintf("std_%s", rand_factor)))[min_index]
min_list = c(rand_factor, min_name, min_std)
}
}
print(min_list)
rownames(get(sprintf("std_%s", rand_factor)))
rand_intercepts = list()
for (rand_factor in rand_factors) {
if (rand_factor == min_list[1]) {
x = rownames(get(sprintf("std_%s", rand_factor)))[-1]
y = x[x != min_list[2]]
if (length(y) == 0) {
rand_intercept = sprintf("(1 |%s)", rand_factor)
}
else {
z = paste(y, collapse = " + ")
rand_intercept = sprintf("(1 + %s|%s)", z, rand_factor)
}
rand_intercepts = append(rand_intercepts, rand_intercept)
}
else {
x = rownames(get(sprintf("std_%s", rand_factor)))[-1]
y = x
if (length(y) == 0) {
rand_intercept = sprintf("(1 |%s)", rand_factor)
}
else {
z = paste(y, collapse = " + ")
rand_intercept = sprintf("(1 + %s|%s)", z, rand_factor)
}
rand_intercepts = append(rand_intercepts, rand_intercept)
}
}
rs = paste(rand_intercepts, collapse = " + ")
formula = as.character(model_summary$call)
lme_function = formula[1]
lme_formula = formula[2]
lme_data = formula[3]
dependent_independent = strsplit(lme_formula, "\(")[[1]][1]
new_line = paste(dependent_independent, rs)
message(paste("The old formula is: ", lme_formula))
print(variances)
message(paste("The old formula is: ", lme_formula))
message(sprintf("removed %s in %s", min_list[2], min_list[1]))
message(paste("The new formula is: ", new_line))
i = i + 1
print(min_list)
if (min_list[2][[1]] %>% is.null) {
if (beeping) {
beep(5)
}
return(ps)
}
old_model_name = deparse(substitute(model))
new_model_name = paste(old_model_name, as.character(i), sep = "_")
str_formula = sprintf("%s = lmer(%s, data=%s)", new_model_name,
new_line, lme_data)
print(str_formula)
eval(parse(text = str_formula))
appender = sprintf("ps = append(ps, %s)", new_model_name)
eval(parse(text = appender))
print(appender)
print(ps)
recall = sprintf("step(%s, i, beeping, ps=ps)",
new_model_name)
eval(parse(text = recall))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.