Heckman procedure on a complex survey data in R Any idea on how to run the Heckman correction on a complex survey data in R?
I've tried doing it manually, but no success so far...
For the first stage I ran the svyglm() function from the survey package which works well and I was able to estimate the probit model . However, for the second stage I'm having trouble including the predicted inverse Mills ratio (λ) in the svyglm() function.

are you looking for something like the two-step approach described here?
https://m-clark.github.io/models-by-example/heckman-selection.html#two-step-approach
library(survey)

set.seed(123456)

N = 10000
educ = sample(1:16, N, replace = TRUE)
age  = sample(18:64, N, replace = TRUE)
wgt = sample( 1:10, N, replace = TRUE)

covmat = matrix(c(.46^2, .25*.46, .25*.46, 1), ncol = 2)
errors = mvtnorm::rmvnorm(N, sigma = covmat)
z = rnorm(N)
e = errors[, 1]
v = errors[, 2]

wearnl = 4.49 + .08 * educ + .012 * age + e

d_star = -1.5 + 0.15 * educ + 0.01 * age + 0.15 * z + v

observed_index  = d_star > 0

d_df = data.frame(wearnl, educ, age, z, observed_index)

d_surv <- svydesign( ~ 1 , data = d_df , weights = ~ wgt )


probit = svyglm(observed_index ~ educ + age + z,
             design   = d_surv,
             family = binomial(link = 'probit'))

probit_lp = predict(probit)

d_surv <- update( d_surv , mills0 = dnorm(probit_lp)/pnorm(probit_lp))

lm_select <- 
    svyglm(
        wearnl ~ educ + age + mills0, 
        design = subset( d_surv , observed_index ) 
    )

summary(lm_select)

