Propensity score matching in cross-sectional survey data and different size of subgroups Let's say that I have a mutation in 1% of a population in a cross-sectional survey:
library(NHANES)
library(survey)
library(MatchIt)

# Generate survey dataset including sample weights and strata
myvars <- c("ID","Gender", "Age", "BMI", "Race1", "Diabetes","SDMVPSU", "SDMVSTRA", "WTMEC2YR")
nhanes <- as.data.frame(subset(NHANESraw, SurveyYr == "2009_10"))
nhanes <- na.omit(nhanes[myvars])

# Including presence of "mutation" in 1% of the entire population. Yes = 1; No = 0
set.seed(999)
nhanes$mutation <- sample(c(0, 1), replace = T, size = nrow(nhanes), prob = c(0.99, 0.01))

Now I would like to examine whether there is a difference in BMI mean according to the presence of the mutation. 
As expected, there is a huge between-group difference in size according to the presence or absence of the mutation:
mutation0 = 9300; mutation1 = 109.
In order to achieve my task I attempted different approaches:
1. Crude unweighted GLM
# crude unweighted GLM, which showed a difference in mean BMI of 1.47236 between groups (P = 0.0482)

summary(glm(BMI ~ mutation, data = nhanes))

2. Crude weighted GLM
# crude weighted GLM, which showed a difference in mean BMI of 2.4454 between groups (P = 0.0579)
# Using sarvey procedure accounting for sample weight
nhanesDesign <- svydesign(id      = ~SDMVPSU,
                          strata  = ~SDMVSTRA,
                          weights = ~WTMEC2YR,
                          nest    = TRUE,
                          data    = nhanes)
summary(svyglm(BMI ~ mutation, 
              data   = nhanes,
              design = nhanesDesign))

3. ADJ unweighted GLM
# ADJ unweighted GLM, which showed a difference in mean BMI of 0.97542 between groups (P =  0.126)
summary(glm(BMI ~ mutation + Age + Gender + Race1 + Diabetes, data = nhanes))

4. ADJ weighted GLM
# ADJ weighted GLM, which showed a difference in mean BMI of 1.831912 between groups (P =  0.0620)
# Using sarvey procedure accounting for sample weight
    summary(svyglm(BMI ~ mutation + Age + Gender + Race1 + Diabetes,
                   data   = nhanes,
                   design = nhanesDesign))

I feel that there is something that I miss. I guess that this huge difference in size between groups can lead to misleading results. 
For this reason I have approched Propensity Score Matching (excluding BMI variable) in order to provide more reliable results:
5. Propensity Score Matching
# Choose subjects according to propensity score
match_model <- matchit(mutation ~ Age + Gender + Race1 + Diabetes, 
                       data = nhanes, 
                       method = "nearest")
match_data <- match.data(match_model)
colnames(match_data)[12] <- "propensity_weights"


# ADJ unweighted GLM from Propensity
summary(glm(BMI ~ mutation + Age + Gender + Race1 + Diabetes, data = match_data))



# Survey Design from Propensity matching
nhanesPropensity <- svydesign(id      = ~SDMVPSU,
                              strata  = ~SDMVSTRA,
                              weights = ~WTMEC2YR,
                              nest    = TRUE,
                              data    = match_data)
summary(svyglm(BMI ~ mutation + Age + Gender + Race1 + Diabetes,
               data   = match_data,
               design = nhanesPropensity))

I have some questions: 


*

*Which one of the above listed methods is the most reliable? In
other words, what is the right approach to data analyses
(specifically, for the calculation of mean BMI) when there is a
large difference in size between groups, like in this case? 

*is propensity score matching an appropriate approach to overcome
between-groups differences in this kind of population
(cross-sectional survey)?

*is the use of sample weights appropriate in propensity matching, as proposed above?

 A: First, you would expect that the mutation and no mutation groups would have the same BMI since the presence of the mutation is completely independent from BMI and all other variables in the data. The fact that you got a significant difference using the first comparison is a fluke; if you use a different seed, (e.g., 9999), you get a nonsignificant difference. No matter what method you use, your estimate of the effect of mutation should be unbiased for 0. That is, whether you consider the unweighted population, the survey-weighted population, or a propensity score-matched population, the difference in group means for BMI is equal to zero, and no amount of conditioning or adjusting will make any of the estimators fail to converge to zero in expectation. The estimator with the lowest variance (i.e., the most precision) will be the covariate-adjusted unweighted comparison because adding covariates predictive of the outcome increases precision and including survey weights decreases precision.
If your question is more broadly about estimating an exposure effect in complex survey data, then we can discuss the differences among the approaches. First, it doesn't matter that your groups are of unequal size. That affects the precision of any estimator, but not the choice of which estimator you should use.
Second, you only need to incorporate the survey design into the analysis if the quantity of interest varies across sampling units or across the variables used to compute the survey weights. If not (i.e., if the exposure effect was equal for all members of the population), you don't need to incorporate the sampling design. Failing to incorporate the survey design means that your estimates generalize to the wrong population, but if the quantity estimated doesn't vary across populations, it doesn't matter whether you attempt to generalize to the unweighted or correctly weighted population.
Third, propensity score analysis with complex survey designs is not well understood, but there has been some literature on it. See DuGoff, Schuler, and Stuart (2014) for a review and current best practices. The choices of how to proceed depend on the goals of the analysis.
Finally, to answer an implied question, it's rarely the case that propensity score methods outperform well-specified regression models. Propensity score methods severely decrease the precision of your estimate, thereby increasing the typical error of an estimate even if the bias is reduced. Unless you have a good reason to prefer propensity score-based approaches (and the one you provided, to generally improved robustness, is not one), you should probably stick with regression approaches to estimate exposure effects, though you should use more sophisticated regression methods like the parametric g-formula over just throwing covariates in a regression.

DuGoff, E. H., Schuler, M., & Stuart, E. A. (2014). Generalizing Observational Study Results: Applying Propensity Score Methods to Complex Surveys. Health Services Research, 49(1), 284–303. https://doi.org/10.1111/1475-6773.12090
