How to obtain standardized betas from a svyglm() regression

Question

I have a svyglm() weighted linear regression model and would like to obtain the standardized betas of the regression coefficients. I have tried to get them with the lm.beta() function, however when I do this the p-values of the coefficients change. I have also tried with the effectsize:standardize_parameters() function, but this returns the unstandardized coefficients (the same ones as in the original svyglm() model). How can I obtain the standardized coefficients of the svyglm() model with the correct p-values? Any help would be greatly appreciated. Here is the code for the model:

library(tidyverse)
library(lm.beta)
library(survey)
library(sjstats)
dat <- read.csv("https://raw.githubusercontent.com/LucasTremlett/questions/master/questiondata.csv")
dat.weighted<- svydesign(ids = ~1, data = dat, weights = dat$weight)
model.weighted<- svyglm(DV~IV1+IV2+IV3, design=dat.weighted)
summary(model.weighted)

Here is the code for the two methods that have not worked for me

##Withlmbeta()
model.weighted.lm.beta <- lm.beta(model.weighted)
summary(model.weighted.lm.beta)

##Witheffectsize:standardize_parameters()
effectsize::standardize_parameters(model.weighted)

Answer 1

I am not sure if lm.beta() can be used for your svyglm object - it's a svy regression and also of glm-type, not lm-type.

If the coefficients returned by lm.beta() are correct (or you can get those form elsewhere), getting standard deviations of the 'standardized beta' seems straightforward.

A common calculation of the 'standardized beta' suggests that it is just multiplied by the ratio of standard deviations of x and y. Then, the standardized beta, a random variable, is also given as the product of this ratio and the estimated s.e. of beta.

# standardized coefficients
st_mean <-  model.weighted.lm.beta$standardized.coefficients 

# standardized coefficients' s.d.
st_sd <- estimates_sd * standadizing_factors

# standardized coefficients tstat
st_tstat <- st_mean/st_sd

result <- data.frame(tibble(st_mean, st_sd, st_tstat))
row.names(result) <- names(model.weighted$coefficients)
result %>% round(5)

Answer 2

Here's a model

data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
model <- svyglm(api00~mobility, design=dstrat)

And we want the estimated population standard deviations

vars<-svyvar(~api00+mobility, design=dstrat)
sd.y <- sqrt(vars[1,1])
sd.x <- sqrt(vars[2,2])

And now the standardised beta

unstd.beta <- coef(model)[2]
std.beta <- unstd.beta*sd.x/sd.y