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