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.

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().
