anova.svyglm.Rd
A method for the anova
function, for use on
svyglm
and svycoxph
objects. With a single model argument it produces a sequential anova table, with two arguments it compares the two models.
# S3 method for svyglm
anova(object, object2 = NULL, test = c("F", "Chisq"),
method = c("LRT", "Wald"), tolerance = 1e-05, ..., force = FALSE)
# S3 method for svycoxph
anova(object, object2=NULL,test=c("F","Chisq"),
method=c("LRT","Wald"),tolerance=1e-5,...,force=FALSE)
# S3 method for svyglm
AIC(object,...,k=2, null_has_intercept=TRUE)
# S3 method for svyglm
BIC(object,...,maximal)
# S3 method for svyglm
extractAIC(fit,scale,k=2,..., null_has_intercept=TRUE)
# S3 method for svrepglm
extractAIC(fit,scale,k=2,..., null_has_intercept=TRUE)
Use (linear combination of) F or chi-squared distributions for p-values. F is usually preferable.
Use weighted deviance difference (LRT) or Wald tests to compare models
For models that are not symbolically nested, the tolerance for deciding that a term is common to the models.
For AIC
and BIC
, optionally more svyglm
objects
not used
Does the null model for AIC have an intercept or not?
Force the tests to be done by explicit projection even if the models are symbolically nested (eg, for debugging)
A svyglm
model that object
(and ... if supplied) are nested in.
Multiplier for effective df in AIC. Usually 2. There is no choice of k
that will give BIC
The reference distribution for the LRT depends on the misspecification effects for the parameters being tested (Rao and Scott, 1984). If the models are symbolically nested, so that the relevant parameters can be identified just by manipulating the model formulas, anova
is equivalent to regTermTest
.If the models are nested but not symbolically nested, more computation using the design matrices is needed to determine the projection matrix on to the parameters being tested. In the examples below, model1
and model2
are symbolically nested in model0
because model0
can be obtained just by deleting terms from the formulas. On the other hand, model2
is nested in model1
but not symbolically nested: knowing that the model is nested requires knowing what design matrix columns are produced by stype
and as.numeric(stype)
. Other typical examples of models that are nested but not symbolically nested are linear and spline models for a continuous covariate, or models with categorical versions of a variable at different resolutions (eg, smoking yes/no or smoking never/former/current).
A saddlepoint approximation is used for the LRT with numerator df greater than 1.
AIC
is defined using the Rao-Scott approximation to the weighted
loglikelihood (Lumley and Scott, 2015). It replaces the usual penalty term p, which is the null expectation of the log likelihood ratio, by the trace of the generalised design effect matrix, which is the expectation under complex sampling. For computational reasons everything is scaled so the weights sum to the sample size.
BIC
is a BIC for the (approximate) multivariate Gaussian models
on regression coefficients from the maximal model implied by each
submodel (ie, the models that say some coefficients in the maximal model
are zero) (Lumley and Scott, 2015). It corresponds to comparing the models with a Wald test and replacing the sample size in the penalty by an effective sample size.
For computational reasons, the models must not only be nested, the names of the coefficients must match.
extractAIC
for a model with a Gaussian link uses the actual AIC based on maximum likelihood estimation of the variance parameter as well as the regression parameters.
Object of class seqanova.svyglm
if one model is given, otherwise of class regTermTest
or regTermTestLRT
At the moment, AIC
works only for models including an intercept.
Rao, JNK, Scott, AJ (1984) "On Chi-squared Tests For Multiway Contingency Tables with Proportions Estimated From Survey Data" Annals of Statistics 12:46-60.
Lumley, T., & Scott, A. (2014). "Tests for Regression Models Fitted to Survey Data". Australian and New Zealand Journal of Statistics, 56 (1), 1-14.
Lumley T, Scott AJ (2015) "AIC and BIC for modelling with complex survey data" J Surv Stat Methodol 3 (1): 1-18.
data(api)
dclus2<-svydesign(id=~dnum+snum, weights=~pw, data=apiclus2)
model0<-svyglm(I(sch.wide=="Yes")~ell+meals+mobility, design=dclus2, family=quasibinomial())
model1<-svyglm(I(sch.wide=="Yes")~ell+meals+mobility+as.numeric(stype),
design=dclus2, family=quasibinomial())
model2<-svyglm(I(sch.wide=="Yes")~ell+meals+mobility+stype, design=dclus2, family=quasibinomial())
anova(model2)
#> Error in .svycheck(design): object 'dclus2' not found
anova(model0,model2)
#> Working (Rao-Scott+F) LRT for stype
#> in svyglm(formula = I(sch.wide == "Yes") ~ ell + meals + mobility +
#> stype, design = dclus2, family = quasibinomial())
#> Working 2logLR = 21.52228 p= 0.0013407
#> (scale factors: 1.7 0.3 ); denominator df= 34
anova(model1, model2)
#> Working (Rao-Scott+F) LRT for stype - as.numeric(stype)
#> in svyglm(formula = I(sch.wide == "Yes") ~ ell + meals + mobility +
#> stype, design = dclus2, family = quasibinomial())
#> Working 2logLR = 25.10744 p= 1.816e-05
#> df=1; denominator df= 34
anova(model1, model2, method="Wald")
#> Wald test for stype - as.numeric(stype)
#> in svyglm(formula = I(sch.wide == "Yes") ~ ell + meals + mobility +
#> stype, design = dclus2, family = quasibinomial())
#> F = 17.73552 on 1 and 34 df: p= 0.00017598
AIC(model0,model1, model2)
#> eff.p AIC deltabar
#> [1,] 3.538287 142.11041 1.179429
#> [2,] 6.248895 119.31208 1.562224
#> [3,] 7.202561 97.03351 1.440512
BIC(model0, model2,maximal=model2)
#> p BIC neff
#> [1,] 4 120.8183 72.43313
#> [2,] 6 111.6461 NaN
## from ?twophase
data(nwtco)
dcchs<-twophase(id=list(~seqno,~seqno), strata=list(NULL,~rel),
subset=~I(in.subcohort | rel), data=nwtco)
a<-svycoxph(Surv(edrel,rel)~factor(stage)+factor(histol)+I(age/12), design=dcchs)
b<-update(a, .~.-I(age/12)+poly(age,3))
## not symbolically nested models
anova(a,b)
#> Working (Rao-Scott+F) LRT for poly(age, 3) - I(age/12)
#> in svycoxph(formula = Surv(edrel, rel) ~ factor(stage) + factor(histol) +
#> poly(age, 3), design = dcchs)
#> Working 2logLR = 9.743054 p= 0.0086592
#> (scale factors: 1.2 0.82 ); denominator df= 1146