I have two periods of panel data and I am trying to see if coefficients change over time. I would like to control for individual-level heterogeneity, but not sure how to do it with a complex design object and svyglm function:

# Libraries
library(survey)

# complex survey design
des <- svydesign(
  ~secu , # stratum half-sample code
  strata = ~stratum , # std error stratum
  weights = ~wtresp , # weights
  nest = TRUE ,
  data = long
)

#
summary(svyglm(
  retired ~ time*(female + college + factor(race) + age + Rshlt + Sshlt + log(iearn + 1) + factor(month) + factor(year)),
  design = des
))

I have tried using plm() with weights, but not sure if it is an appropriate method for complex survey samples.

Thank you!