- dfP <- sum(phy$edge.length)*N / sum(diag(vcv.phylo(phy)))
- obj <- list(call = geemod$call,
+ ## <FIXME>
+ ## maybe need to refine below in case of non-Brownian corStruct
+ if (!missing(corStruct)) phy <- attr(corStruct, "tree")
+ dfP <- sum(phy$edge.length)*N / sum(diag(R))
+ ## </FIXME>
+
+ ## compute QIC:
+ Y <- geemod$y
+ MU <- geemod$fitted.values
+ Qlik <- switch(fname,
+ "gaussian" = -sum((Y - MU)^2)/2,
+ "binomial" = sum(Y*log(MU/(1 - MU)) + log(1 - MU)),
+ "poisson" = sum(Y*log(MU) - MU),
+ "Gamma" = sum(Y/MU + log(MU)),
+ "inverse.gaussian" = sum(-Y/(2*MU^2) + 1/MU))
+ Ai <- do.call("gee", list(formula, id, data = data, family = family,
+ corstr = "independence", scale.fix = scale.fix,
+ scale.value = scale.value))$naive.variance
+ QIC <- -2*Qlik + 2*sum(diag(solve(Ai) %*% W))
+
+ obj <- list(call = match.call(),