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.

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.
