-
Notifications
You must be signed in to change notification settings - Fork 2
/
new_lmer_AIC_tables3.r
149 lines (145 loc) · 4.83 KB
/
new_lmer_AIC_tables3.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
aicW <- function(model.list, finite = TRUE, null.model = NULL, N, order = TRUE){
if(hasArg(N)) { N <- N }
else N <- NULL
if(is.null(null.model)) {
null.model <- model.list[[length(model.list)]]
}
LLlist <- vector()
for(index in 1:length(model.list)){
if(class(model.list[[index]])[1] == "mer") {
LLlist <- c(LLlist, as.numeric(logLik(model.list[[index]])))
}
else {LLlist <- c(LLlist, as.numeric(logLik(model.list[[index]])))}
}
aiclist <- vector()
for(index in 1:length(model.list)){
aicname <- "AIC"
if(finite==T) {
aiclist <- c(aiclist, AIC.finite.lm(model.list[[index]], N))
aicname <- "AICc"}
else {aiclist <- c(aiclist, AIC.lm(model.list[[index]]))}
}
biclist <- vector()
for(index in 1:length(model.list)){
if(class(model.list[[index]])[1] == "mer") {
biclist <- c(biclist, BIC.lm(model.list[[index]], N))
}
else {biclist <- c(biclist, BIC.lm(model.list[[index]], N))}
}
dflist <- vector()
for(index in 1:length(model.list)){
if(class(model.list[[index]])[1] == "mer") {
dflist <- c(dflist, as.numeric(attributes(logLik(model.list[[index]]))$df))
}
else {dflist <- c(dflist, length(coef(model.list[[index]])))}
}
delist <- vector()
for(index in 1:length(model.list)){
if(class(model.list[[index]])[1] == "mer") {
delist <- c(delist, lmer.de(model.list[[index]], null.model, N))
}
else {delist <- c(delist, lmer.de(model.list[[index]], null.model, N))}
}
modform <- vector()
if(is.null(names(model.list))){
for(index in 1:length(model.list)){
modform <- c(modform, formula(model.list[[index]]))
}
}
if(!is.null(names(model.list))){
models <- names(model.list)
delta.i <- aiclist - min(aiclist)
aicw <- exp(delta.i * -0.5)/sum(exp(delta.i * -0.5))
delta.BIC <- biclist - min(biclist)
bicw <- exp(delta.BIC * -0.5)/sum(exp(delta.BIC * -0.5))
return.matrix <- matrix(c(round(LLlist, 3), dflist, round(aiclist, 3), round(delta.i, 3), round(aicw, 4), round(biclist,3), round(delta.BIC, 3), round(bicw,4), round(delist, 2)), ncol = 9)
colnames(return.matrix) <- c("logLik", "df", aicname, paste("d", aicname, sep = ""), "wAICc", "BIC", "dBIC", "wBIC", "%DE")
# rownames(return.matrix) <- sapply(strsplit(models, NULL), function(x) paste(x, collapse = " + "))
rownames(return.matrix) <- models
}
else {
delta.i <- aiclist - min(aiclist)
aicw <- exp(delta.i * -0.5)/sum(exp(delta.i * -0.5))
delta.BIC <- biclist - min(biclist)
bicw <- exp(delta.BIC * -0.5)/sum(exp(delta.BIC * -0.5))
return.matrix <- matrix(c(round(LLlist, 3), dflist, round(aiclist, 3), round(delta.i, 3), round(aicw, 4), round(biclist,3), round(delta.BIC, 3), round(bicw,4), round(delist, 2)), ncol = 9)
# return.matrix <- matrix(c(round(LLlist, 3), dflist, round(aiclist, 3), round(delta.i, 3), round(aicw, 4), round(delta.BIC, 3), round(delist, 2)), ncol = 7)
colnames(return.matrix) <- c("logLik", "df", aicname, paste("d", aicname, sep = ""), "wAICc", "BIC", "dBIC", "wBIC", "%DE")
# colnames(return.matrix) <- c("logLik", "df", aicname, paste("d", aicname, sep = ""), "weight", "dBIC", "%DE")
rownames(return.matrix) <- modform
}
if(order) {
return.matrix <- return.matrix[order(return.matrix[, 4]), ]
}
return(return.matrix)
}
AIC.finite.lm <- function(model, N){
LL = as.numeric(logLik(model))
p = as.numeric(attributes(logLik(model))$df)
if(class(model)[1] == "mer"){
if(is.null(N)){
N <- model@A@Dim[2]
}
else N <- N
}
else {
if(is.null(N)){
N <- length(model$data[, 1])
}
else N <- N
}
AICc <- AIC(logLik(model)) + (2 * p * (p + 1))/(N - p - 1)
return(AICc)
}
AIC.lm <- function(object) {
ret <- AIC(logLik(object))
return(ret)
}
BIC.lm <- function(model, N) {
LL = as.numeric(logLik(model))
p = as.numeric(attributes(logLik(model))$df)
if(class(model)[1] == "mer"){
if(is.null(N)){
N <- model@A@Dim[2]
}
else N <- N
}
else {
if(is.null(N)){
N <- length(model$data[, 1])
}
else N <- N
}
BIC <- (-2 * LL) + (p * log(N))
return(BIC)
}
lmer.de <- function(model, null.model, N) {
if(class(model)[1] == "mer"){
if(is.null(N)){
N <- model@A@Dim[2]
}
else N <- N
}
else {
if(is.null(N)){
N <- length(model$data[, 1])
}
else N <- N
}
if(class(null.model)[1] == "mer"){
null.dev <- as.vector(null.model@deviance["ML"])
model.lr <- as.vector(null.model@deviance["ML"]) - as.vector(model@deviance["ML"])
}
else {
if(class(model)[1] == "mer"){
null.dev <- deviance(null.model)
model.lr <- (deviance(null.model)) - as.vector(model@deviance["ML"])
}
else {
null.dev <- deviance(null.model)
model.lr <- deviance(null.model) - deviance(model)
}
}
de <- 100 * (model.lr/null.dev)
return(de)
}