1
2
3logLik.svyglm<-function(object,...){
4   warning("svyglm not fitted by maximum likelihood.")
5   object$deviance/-2
6}
7
8AIC.svyglm<-function(object,...,k=2,null_has_intercept=TRUE){
9	if (length(list(...))){
10		do.call(rbind,lapply(list(object,...),extractAIC,k=k,null_has_intercept=null_has_intercept))
11    } else {
12	   extractAIC(object,k=k,null_has_intercept=null_has_intercept)
13    }
14}
15extractAIC.svyglm<-function(fit,scale,k=2,...,null_has_intercept=TRUE){
16    if(is.svylm(fit)) return(extractAIC_svylm(fit,...))
17    if (length(attr(terms(fit),"factors"))){
18        ftest<-delete.response(formula(fit))
19        if (!null_has_intercept)
20            ftest<-update(ftest,.~.+`1`)
21        r<-regTermTest(fit, ftest, method="LRT")
22        deltabar<-mean(r$lambda)
23	} else {
24	    r<-list(lambda=0)
25	    deltabar<-NaN
26	}
27	d<-fit$deviance
28	c(eff.p=sum(r$lambda), AIC=d+k*sum(r$lambda),deltabar=deltabar)
29}
30
31extractAIC.svrepglm<-extractAIC.svyglm
32
33BIC.svyglm<-function(object,...,maximal){
34	if (length(list(...))){
35		do.call(rbind,lapply(list(object,...),dBIC,modelM=maximal))
36    } else {
37	   dBIC(object,modelM=maximal)
38    }
39
40}
41
42dBIC<-function(modela, modelM){
43	pm<-modela$rank
44	pM<-modelM$rank
45
46	if (any(!(names(coef(modela))%in% names(coef(modelM))))){
47		stop("coefficients in model but not in maximal model")
48		}
49	index<-!(names(coef(modelM))%in% names(coef(modela)))
50	n<-1+modela$df.null
51	if(any(index)){
52		wald<-coef(modelM)[index]%*%solve(vcov(modelM)[index,index],coef(modelM)[index])
53		detDelta<-det(solve(modelM$naive.cov[index,index,drop=FALSE],modelM$cov.unscaled[index,index,drop=FALSE]))
54		dbar<-detDelta^(1/(pM-pm))
55		nstar<-n/dbar
56	}else {
57		wald<-0
58		detDelta<-1
59		dbar<-1
60		nstar=NaN
61		}
62	c(p=pm, BIC=wald+pm*log(n)+log(detDelta)+deviance(modelM),neff=nstar)
63    }
64
65extractAIC.svrepcoxph<-function (fit, scale, k = 2, ...) .NotYetImplemented()
66extractAIC.svycoxph<-function (fit, scale, k = 2, ...)
67{
68    Delta<-solve(fit$inv.info, fit$var)
69    deltabar <- mean(diag(Delta))
70    d <- -2*fit$ll[1]
71    c(eff.p = sum(diag(Delta)), AIC = d + k * sum(diag(Delta)), deltabar = deltabar)
72}
73AIC.svycoxph<-function (object, ..., k = 2)
74{
75    if (length(list(...))) {
76        do.call(rbind, lapply(list(object, ...), extractAIC,
77            k = k))
78    }
79    else {
80        extractAIC(object, k = k)
81    }
82}
83
84
85##  special-case the Gaussian
86## dAIC=-2\max_{\beta,\sigma^2}\log\ell +2p= n\log\mathrm{RSS}-n\log n +n +n\log 2\pi +2p\bar\delta
87
88is.svylm<-function(it) {inherits(it,"svyglm") && isTRUE(all.equal(stats::family(it),gaussian()))}
89
90extractAIC_svylm<-function(fit,scale,k=2,...,null_has_intercept=TRUE){
91    sfit <- summary(fit)
92    n<-sfit$df.null
93    Nhat<-sum(w<-fit$prior.weights)
94    sigma2hat<-coef(sfit$dispersion)
95
96    minus2ellhat<- Nhat*log(sigma2hat) +Nhat +Nhat*log(2*pi)
97
98    V0<-fit$naive.cov
99    V<-vcov(fit)
100    if (null_has_intercept){
101        V0<-V0[-1,-1,drop=FALSE]
102        V<-V[-1,-1,drop=FALSE]
103    }
104
105    ## for mu
106    Delta_mu<-solve(V0,V)
107
108    ## Now for sigma2
109    y<-fit$y
110    muhat<-fit$linear.predictors
111    Isigma2<-Nhat/(2*sigma2hat^2)
112    Usigma2<- -1/(2*sigma2hat)+ (y-muhat)^2/(2*sigma2hat^2)
113    Hsigma2<- sum(w*Usigma2^2)
114    Deltasigma2<-Hsigma2/Isigma2
115
116    ## combine
117    deltabar<-mean(c(diag(Delta_mu), Deltasigma2))
118    eff.p<-sum(diag(Delta_mu))+Deltasigma2
119    aic <- minus2ellhat + k*eff.p
120
121    c(eff.p=eff.p, AIC=aic,deltabar=deltabar)/sigma2hat
122}
123
124
125