# Test support for different classes of models

require(MuMIn)
.checkPkg <- function(package) length(.find.package(package, quiet=TRUE)) != 0

# TEST gls --------------------------------------------------------------------------------
library(nlme)

fm1Dial.gls <- gls(rate ~(pressure + I(pressure^2) + I(pressure^3))*QB,
      Dialyzer, method="ML")

varying <- list(
	correlation = alist(
		AR1_0.771=corAR1(0.771, form = ~ 1 | Subject), 	AR1=corAR1(), NULL),
	weights = alist(vp.press=varPower(form = ~ pressure), NULL)
	)

dd <- dredge(fm1Dial.gls, m.max = 2, m.min = 1, fixed=~pressure, varying = varying, extra = "R^2")

models <- get.models(dd, 1:4)


ma <- model.avg(models, revised=T)
ms <- model.sel(models)
print(ms, abbr = F)
print(ms, abbr = T)
summary(ma)
predict(ma)[1:10]

predict(ma, newdata = Dialyzer[1:5, ])

detach(package:nlme); rm(list=ls())

# TEST glmmML --------------------------------------------------------------------
if(.checkPkg("glmmML")) {
library(glmmML)

set.seed(100)
dat <- data.frame(y = rbinom(100, prob = rep(runif(20), rep(5, 20)), size = 1),
	x = rnorm(100), x2 = rnorm(100), id = factor(rep(1:20, rep(5, 20))))

fm1 <- glmmML(y ~ x*x2, data = dat, cluster = id)
dd <- dredge(fm1)
# mod <- get.models(dd, subset = delta <= 4)
summary(ma <- model.avg(dd, subset = delta <= 4))
# vcov(ma)
coefTable(ma)

detach(package:glmmML); rm(list=ls())
}

# TEST lm ---------------------------------------------------------------------------------
if(.checkPkg("nlme")) {
data(Orthodont, package = "nlme")

fm1 <- lm(distance ~ Sex*age + age*Sex, data = Orthodont)

dispersion <- function(object) {
	wts <- weights(object); if(is.null(wts)) wts <- 1
	sum((wts * resid(object, type = "working")^2)[wts > 0]) / df.residual(object)
}

dd <- dredge(fm1, extra = alist(dispersion))
gm <- get.models(dd, 1:4)
ma <- model.avg(gm, revised=F)

vcov(ma)
summary(ma)
confint(ma)

predict(ma)
predict(ma, se.fit=TRUE)
predict(ma, data.frame(Sex="Male", age=8:12))

rm(list=ls())
}

# TEST glm --------------------------------------------------------------------------------
data(Cement, package = "MuMIn")

nseq <- function(x, len=length(x)) seq(min(x, na.rm=TRUE), max(x, na.rm=TRUE),
	length=len)

fm1 <- glm(y ~ (X1+X2+X3)^2, data = Cement)

dd <- dredge(fm1)
gm <- get.models(dd, 1:10)

summary(ma <- model.avg(gm))
#summary(ma1<-model.avg(dd, 1:10))
#summary(ma2<-model.avg(dd[1:10]))


#vcov(fm1)
vcov(ma)

summary(ma1  <- model.avg(dd[1:10]))
summary(ma2 <- model.avg(model.sel(dd[1:10], rank = "AICc")))

all.equal(ma$avg.model, ma1$avg.mode)

predict(ma) == predict(ma, Cement)
predict(ma, se.fit=T)
predict(ma, lapply(Cement, nseq))


rm(list=ls())
# TEST rlm --------------------------------------------------------------------------------

if (.checkPkg("MASS")) {
library(MASS)
data(Cement, package = "MuMIn")

nseq <- function(x, len=length(x)) seq(min(x, na.rm=TRUE), max(x, na.rm=TRUE),
	length=len)

fm1 <- rlm(y ~X1+X2*X3+X4, data = Cement)
dd <- dredge(fm1, trace=T)
gm <- get.models(dd, 1:10)
ma <- model.avg(gm)
stopifnot(all(predict(ma) == predict(ma, Cement)))
predict(ma, lapply(Cement, nseq, len=30), se.fit=TRUE)
vcov(ma)

rm(list=ls()); #detach(package:MASS)

# TEST multinom --------------------------------------------------------------------------------
if (.checkPkg("nnet")) {
library(nnet); library(MASS)

# Trimmed-down model from example(birthwt)
data(birthwt)

bwt <- with(birthwt, data.frame(
		low = low,
		race = factor(race, labels = c("white", "black", "other")),
		ptd = factor(ptl > 0),
		smoke = (smoke > 0)
		))

options(contrasts = c("contr.treatment", "contr.poly"))
bwt.mu <- multinom(low ~ ., data = bwt)
dd <- dredge(bwt.mu, trace=T)

summary(model.avg(dd[1:5]))
gm <- get.models(dd, 1:5)
ma <- model.avg(gm)

summary(ma)

# predict(ma) // Cannot average factors!

rm(list=ls()); detach(package:nnet)
}}


# TEST gam --------------------------------------------------------------------------------
if (.checkPkg("mgcv")) {

suppressPackageStartupMessages(library(mgcv))
RNGkind("Mersenne")
set.seed(0) ## simulate some data...
dat <- gamSim(1, n = 400, dist = "binary", scale = 2)
#gam1 <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat)

ops <- options(warn = -1)

gam1 <- gam(y ~ s(x0) + s(x1) + s(x2) +  s(x3) + (x1+x2+x3)^2,
	data = dat, method = "GCV.Cp", family = binomial)

dd <- dredge(gam1, subset=!`s(x0)` & (!`s(x1)` | !x1) & (!`s(x2)` |
	!x2) & (!`s(x3)` | !x3), fixed = "x1")

gm <- get.models(dd, cumsum(weight) <= .95)
ma <- model.avg(gm)

summary(ma)

predict(ma, dat[1:10, ], se.fit=T, type = "link")

predict(ma, dat[1:10, ], se.fit=T, type = "response")
predict(ma, dat[1:10, ], se.fit=T, type = "link", backtransform = TRUE)

options(ops)

rm(list=ls()); detach(package:mgcv)
}

# TEST spautolm ---------------------------------------------------------------------------

if (.checkPkg("spdep"))
if(!is.null(tryCatch(suppressPackageStartupMessages(library(spdep)), error = function(e) NULL))) {

suppressMessages(example(NY_data, echo = FALSE))

# method argument changed in spdep 0.5.53 from "full" to "eigen"
method1 <- if(packageVersion("spdep") < '0.5.53') "full" else "eigen"


fm1.spautolm <- spautolm(Z ~ PEXPOSURE * PCTAGE65P + PCTOWNHOME,
 data = nydata, listw = listw_NY, family = "SAR", method = method1, verbose = FALSE)

options(warn=1)
dd <- dredge(fm1.spautolm, m.max=1, fixed = ~PEXPOSURE,
	varying = list(
		family = list("CAR", "SAR"),
		method=list("Matrix_J", method1)
	), trace=FALSE)
options(warn=0)


#dd <- dredge(fm1.spautolm, m.max=3, fixed=~PEXPOSURE)
gm <- get.models(dd, cumsum(weight) <= .99)
ma <- model.avg(gm)
summary(ma)
signif(resid(ma), 5)[1:10]

rm(list=ls())

# TEST spautolm ---------------------------------------------------------------------------
#suppressPackageStartupMessages(library(spdep))
data(oldcol)

fm1.sarlm <- errorsarlm(CRIME ~ INC * HOVAL * OPEN, data = COL.OLD,
 nb2listw(COL.nb, style = "W"), method = "eigen", quiet = TRUE)

dd <- dredge(fm1.sarlm)

gm <- get.models(dd, cumsum(weight) <= .98)
ma <- model.avg(gm)

stopifnot(isTRUE(all.equal(coefTable(ma), coefTable(model.avg(dd, cumsum(weight) <= .98)))))

summary(ma)

predict(ma)[1:10]

rm(list=ls()); detach(package:spdep)

} # library(spdep)

# TEST glm.nb ---------------------------------------------------------------------------
if (.checkPkg("MASS")) {
require(MASS)

quine.nb1 <- glm.nb(Days ~ 0 + Sex/(Age + Eth*Lrn), data = quine)
#quine.nb1 <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine)

ms <- dredge(quine.nb1, marg.ex = "Sex")

models <- get.models(ms )
summary(model.avg(models))

#dredge(quine.nb1, marg.ex = NULL) # OK
#dredge(quine.nb1, marg.ex = NA) # OK
#dredge(quine.nb1, marg.ex = "Sex") # OK
dredge(quine.nb1) # OK
#dredge(quine.nb1, marg.ex = "Sex") # Right, should be the same as above
ma <- model.avg(dredge(quine.nb1, marg.ex = "Sex"), subset = cumsum(weight)<=.9999)

# Cannot predict with this 'averaging'
#pred <- predict(ma, se=TRUE)

#pred <- cbind(pred$fit, pred$fit - (2 * pred$se.fit), pred$fit + (2 * pred$se.fit))
#matplot(pred, type="l")
#matplot(family(quine.nb1)$linkinv(pred), type="l")

rm(list=ls()); #detach(package:MASS)
}

# TEST quasibinomial -----------------------------------------------------------

budworm <- data.frame(ldose = rep(0:5, 2), numdead = c(1, 4, 9, 13, 18, 20, 0,
	2, 6, 10, 12, 16), sex = factor(rep(c("M", "F"), c(6, 6))))
budworm$SF = cbind(numdead = budworm$numdead, numalive = 20 - budworm$numdead)

qbinomial <- function(...) {
	res <- quasibinomial(...)
	res$aic <- binomial(...)$aic
	res
}
budworm.qqlg <- glm(SF ~ sex*ldose + sex*I(ldose^2), data = budworm, family = qbinomial)
#budworm.qlg <- glm(SF ~ sex*ldose + sex*I(ldose^2), data = budworm, family = quasibinomial)
budworm.lg <- glm(SF ~ sex*ldose + sex*I(ldose^2), data = budworm, family = binomial)

dd <- dredge(budworm.lg, rank = "QAIC", chat = summary(budworm.lg)$dispersion)
#dd <- dredge(budworm.lg) # should be the same
mod <- get.models(dd, seq(nrow(dd)))

# Note: this works:
# ma <- model.avg(mod)
# but this will not: ('rank' attribute passed from 'dredge' is lost)
# ma <- model.avg(mod)
# so, need to supply them
ma <- model.avg(mod[1:5], rank="QAICc", rank.args = list(chat = 0.403111))

rm(list=ls())

# TEST polr {MASS} -------------------------------------------------------------
#if (.checkPkg("MASS")) {
#library(MASS)
#library(MuMIn)
#options(contrasts = c("contr.treatment", "contr.poly"))
#house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)

#dd <- dredge(house.plr)
#mod <- get.models(dd, 1:6)
#model.avg(mod)
#}


# TEST coxph -------------------------------------------------------------------

library(survival)

bladder1 <- bladder[bladder$enum < 5, ]

fmcph <- coxph(Surv(stop, event) ~ (rx + size + number) * strata(enum) +  cluster(id), bladder1)

r.squared.coxph <- function(object, ...) {
	logtest <- -2 * (object$loglik[1L] - object$loglik[2L])
	c(rsq = 1 - exp(-logtest/object$n), maxrsq = 1 - exp(2 * object$loglik[1L]/object$n))
}

ms <- dredge(fmcph, fixed=c("cluster(id)", "strata(enum)"), extra = list(R2="r.squared.coxph"))

fits <- get.models(ms, delta < 5)
summary(model.avg(fits))


fmsrvrg <- survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist='weibull',
    scale = 1)
summary(model.avg(dredge(fmsrvrg), delta  < 4))

rm(list=ls())
detach(package:survival)

# END TESTS

###model1 <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp,
###              data = esoph, family = binomial())
###model1 <- glm(cbind(ncases, ncontrols) ~ agegp,
###              data = esoph, family = binomial())
###model2 <- update(model1, contrasts = list(agegp="contr.treatment"))
###model3 <- update(model1, .~. + tobgp, contrasts = list(agegp="contr.treatment"))
###model3$contrasts
###
###models <- list(model1, model2, model3)
###x <- rbindDataFrameList(lapply(sapply(models, getElement, "contrasts"), as.data.frame))
###
###as.list(x)
###
###any(apply(x, 2, function(x) length(na.omit(unique(x)))) > 1)
###
####system.time(for(i in 1:5000) lapply(lapply(lapply(x, unique), na.omit), length))
####system.time(for(i in 1:5000) lapply(x, function(x) length(na.omit(unique(x)))))
###
###summary(model.avg(model1, model2, model3))
###
###
###model1$contrasts
###
###model2$contrasts
###
###
###ms1 <- dredge(model1, rank=Cp)
###
####plot(dredge(model1, rank=ICOMP))
###dredge(model1, rank=BIC, extra="R^2")
###
###model1$contrasts
###summary(model.avg(ms1, delta <=10))
###
###summary(model2 <- glm(case ~ age+parity+education+spontaneous+induced,
###                data=infert, family = binomial()))
###
###dredge(model2, rank=BIC, extra="nobs")

# require(survival)
# model3 <- clogit(case~spontaneous+induced+strata(stratum),data=infert)
# model3 <- clogit(case~spontaneous+strata(stratum),data=infert)

#require(quantreg)
#data(stackloss)
#fm1 <- rq(stack.loss ~ Water.Temp + Acid.Conc. + Air.Flow, tau = .5, data = stackloss)
#ms <- dredge(fm1, varying = list(tau=list(.75, .5, .25)), extra="R^2")
#fm1 <- lm(stack.loss ~ Water.Temp + Acid.Conc. + Air.Flow,  data = stackloss)
#models <- get.models(ms, delta <= 4)
#model.avg(ms, delta <= 4)
#plot(ms)

#
#library(brglm)
#data(lizards)
#lizg <- glm(cbind(grahami, opalinus) ~ (height + diameter +
#    light + time)^2, family = binomial(logit), data=lizards)
#
#glm(cbind(grahami, opalinus) ~ (height + diameter +
#    light + time)^3, family = poisson("log"), data=lizards)
#dd <- dredge(lizg)
#plot(dd, col=6:7)

#brglm(cbind(grahami, opalinus) ~ height + diameter + light + time,
#	family = binomial(logit), data=lizards, method = "brglm.fit")
#
#lizards.brglm <- brglm(cbind(grahami, opalinus) ~ height + diameter +
#                  light + time, family = binomial(logit), data=lizards,
#                  method = "brglm.fit")