Post hoc analysis (interaction, multiple comparision, ect.) of regression fitted by```svycoxph``` Recently I am working with interaction in regression. And I learn a lot from this
The package emmeans perform well in various kinds of regression.
However, I found it doesn't work with Complex Survey Design data.
So my question is

*

*I would like to conduct multiple comparisons of a multilevel categorical variable in the model fitted by svycoxph like this I used to do

model<-coxph(Surv(futime,death)~factor(BMI_group),data=data) # BMI_group is a 3 level factor
cat<-emmeans(model,~ factor(BMI_group),type = "response")
contrast(cat,'revpairwise',adjust = 'tukey')



*I would like to calculate the main/simple effect in the model with interaction. Below is the example in coxph
model<-coxph(Surv(futime,death)~factor(BMI_group)*factor(Sex),data=data) # BMI_group is a 3 level factor, Sex is a binary variable
catcat<-emmeans(model,~ factor(BMI_group)*factor(Sex),type = "response")
contrast(catcat,'revpairwise',adjust = 'tukey')

When using the Complex Survey Design data, the codes are as below
model<-svycoxph(Surv(futime,Death)~factor(BMI_group),design = Svydata)

AND
model<-svycoxph(Surv(futime,Death)~factor(BMI_group)*factor(Sex),design = Svydata)

Unfortunately package emmeans doesn't work
In a word, I would like to know if there is a substitute for emmeans to deal with Complex Survey Design data. Or if I can calculate what I want manually?
I hope my expression is clear enough and any insight is appreciated.
 A: emmeans() almost works on these models...
Consider one of the examples for survey::svycoxph:
data(pbc, package="survival")

pbc$randomized<-with(pbc, !is.na(trt) & trt>0)
biasmodel<-glm(randomized~age*edema,data=pbc,family=binomial)
pbc$randprob<-fitted(biasmodel)
if (is.null(pbc$albumin)) pbc$albumin<-pbc$alb ##pre2.9.0

dpbc<-svydesign(id=~1, prob=~randprob, strata=~edema, data=subset(pbc,randomized))
rpbc<-as.svrepdesign(dpbc)

(model<-svycoxph(Surv(time,status>0)~log(bili)+protime+albumin,design=dpbc))

If you try
emmeans(model, "bili", at = list(bili = 2:4))

You get an error message saying it can't reconstruct the data. But that error message also suggests providing the dataset as an additional argument. In this case, be a bit careful because this example actually uses a subset of the original dataset; but it works:
> emmeans(model, "bili", at = list(bili = 2:4), 
+         data = subset(pbc, randomized))
 bili  emmean   SE  df asymp.LCL asymp.UCL
    2 -0.4309 1.12 Inf     -2.64      1.77
    3 -0.0717 1.12 Inf     -2.27      2.13
    4  0.1832 1.12 Inf     -2.01      2.38

Results are given on the log (not the response) scale. 
Confidence level used: 0.95

There is a slight additional complication if you want to get results on the response scale (back-transformed from the log scale). If you ask for type = "response" it balks because the model formula has Surv(time,status>0) as the response and it tries to undo that transformation too. So use type = "unlink" instead, which undoes the log link but doesn't try to undo the repsonse transformation.
> emmeans(model, "bili", at = list(bili = 2:4), 
+         data = subset(pbc, randomized), type = "unlink")
 bili response    SE  df asymp.LCL asymp.UCL
    2    0.650 0.731 Inf    0.0717      5.89
    3    0.931 1.044 Inf    0.1033      8.39
    4    1.201 1.346 Inf    0.1336     10.79

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

So in summary, just provide the dataset, and be careful if you want results on the response scale.
A: You can calculate what you want manually, although it's safer to take advantage of tools like those provided by the emmeans package to avoid errors.
The coefficient covariance matrix is the key for evaluating hypotheses based on regression coefficients. For this type of model, the coefficient estimates are assumed to have a multivariate normal distribution that includes correlations among the estimates. The standard errors reported in a model summary() are the square roots of the diagonal values of that matrix.
You need to use the off-diagonal elements of that matrix for many tests, even though they are typically hidden from the model summary(). For example, a test on a linear combination of coefficients (e.g., predictions at particular covariate values) would use the coefficient covariances along with the formula for the variance of a weighted sum of correlated variables. Tests on multiple coefficients at once can use the covariance matrix to perform Wald tests. Apply appropriate corrections for multiple comparisons.
Happily, if your model returns a coefficient covariance matrix via vcov() and a list of coefficient estimates via coef(), then the qdrg() function in emmeans can get those directly from the model, setting up a reference grid that lets other emmeans functions do all that work for you. If vcov() doesn't return a value but you have another way to get that matrix, you can specify it in the vcov argument to qdrg().
