Odette Gaston <odette.gaston <at> gmail.com> writes:
>
> Dear all,
>
> I have a problem with accessing class attributes.
> I was unable to solve this
> yet, but someone may know how to solve it.
My best guess at your immediate problem (doing
things by hand) is that you're not using the
whole vector. From your example:
Delta <- c(m1 = 0, m2 = 1.8, m3 = 4.2, m4 = 6.2)
exp(-0.5*Delta)/sum(exp(-0.5*Delta))
m1 m2 m3 m4
0.63529363 0.25829111 0.07779579 0.02861947
In general the dRedging package at
http://www.zbs.bialowieza.pl/users/kamil/r/ can do these
problems (I hate to recommend this package because it
offers the danger of thoughtless convenience,
but if you really know that you want to enumerate
models and do IC-based model averaging it can save a
lot of time). At the moment, though, it doesn't work
with glmmML-based objects (you could ask the author
to extend it).
When I tried stepAIC it didn't really enumerate
all the models for me (that's not its purpose),
so I went through and enumerated by hand. For example;
library(glmmML)
set.seed(1001)
a <- runif(100)
b <- runif(100)
c <- runif(100)
x <- runif(100)
n <- rep(20,100)
cluster <- factor(rep(1:5,20))
linpred <- a+b+c+x-2
y <- rbinom(100,prob=plogis(linpred),size=n)
data <- data.frame(y,a,b,c,x,n)
m <- list()
## full model
m[[1]] <- glmmML(cbind (y, n-y)~ x+a+b+c,
family = binomial, data, cluster)
## 3-term models
m[[2]] <- update(m[[1]],.~.-a) ## xbc
m[[3]] <- update(m[[1]],.~.-b) ## xac
m[[4]] <- update(m[[1]],.~.-c) ## xab
m[[5]] <- update(m[[1]],.~.-x) ## abc
## 2-term models
m[[6]] <- update(m[[2]],.~.-x) ## bc
m[[7]] <- update(m[[2]],.~.-b) ## xc
m[[8]] <- update(m[[2]],.~.-c) ## xb
m[[9]] <- update(m[[3]],.~.-x) ## ac
m[[10]] <- update(m[[3]],.~.-c) ## xa
m[[11]] <- update(m[[4]],.~.-x) ## ab
## 0-term models (intercept)
m[[12]] <- glmmML(cbind (y, n-y)~ 1, family = binomial, data, cluster)
m[[13]] <- update(m[[12]],.~.+a)
m[[14]] <- update(m[[12]],.~.+b)
m[[15]] <- update(m[[12]],.~.+c)
m[[16]] <- update(m[[12]],.~.+x)
## have to define logLik and AIC for glmmML objects
logLik.glmmML <- function(x) {
loglik <- (-x$deviance)/2
attr(loglik,"df") <- length(coef(x))
loglik
}
AIC.glmmML <- function(x) x$aic
library(bbmle)
## now it works (the answers are pretty trivial
## in this made-up case
AICtab(m,sort=TRUE,weights=TRUE,delta=TRUE)
