Predictive marginal means for a generalised linear model, using the method of Korn and Graubard (1999) and matching the results of SUDAAN. The predictive marginal mean for one level of a factor is the probability-weighted average of the fitted values for the model on new data where all the observations are set to that level of the factor but have whatever values of adjustment variables they really have.

svypredmeans(adjustmodel, groupfactor, predictat=NULL)

## Arguments

A generalised linear model fit by svyglm with the adjustment variable but without the factor for which predictive means are wanted

groupfactor

A one-sided formula specifying the factor for which predictive means are wanted. Can use, eg, ~interaction(race,sex) for combining variables. This does not have to be a factor, but it will be modelled linearly if it isn't

predictat

A vector of the values of groupfactor where you want predictions. If groupfactor is a factor, these must be values in the data, but if it is numeric you can interpolate/extrapolate

## Value

An object of class svystat with the predictive marginal means and their covariance matrix.

## Note

It is possible to supply an adjustment model with only an intercept, but the results are then the same as svymean

It makes no sense to have a variable in the adjustment model that is part of the grouping factor, and will give an error message or NA.

svyglm

Worked example using National Health Interview Survey data: https://gist.github.com/tslumley/2e74cd0ac12a671d2724

## Examples

data(nhanes)
nhanes_design <- svydesign(id=~SDMVPSU, strata=~SDMVSTRA, weights=~WTMEC2YR, nest=TRUE,data=nhanes)
agesexmodel<-svyglm(HI_CHOL~agecat+RIAGENDR, design=nhanes_design,family=quasibinomial)
## high cholesterol by race/ethnicity, adjusted for demographic differences
means<-svypredmeans(agesexmodel, ~factor(race))
means
#>       mean     SE
#> 2 0.114596 0.0065
#> 3 0.084718 0.0106
#> 1 0.123081 0.0058
#> 4 0.108770 0.0304
## relative risks compared to non-Hispanic white
svycontrast(means,quote(1/2))
#>          nlcon     SE
#> contrast 1.074 0.0722
svycontrast(means,quote(3/2))
#>            nlcon     SE
#> contrast 0.73928 0.0923

data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
demog_model <- svyglm(api00~mobility+ell+hsg+meals, design=dstrat)
svypredmeans(demog_model,~enroll, predictat=c(100,300,1000,3000))
#>        mean      SE
#> 100  699.19 10.3203
#> 300  684.29  9.5242
#> 1000 632.13  9.3377
#> 3000 483.11 22.8229