Chapter 3 Analysis of continuous variables
The statistical analysis of surveys has constantly evolved, driven by the development of new methodologies and approaches that seek to improve the precision and representativeness of the results. In general, these advances emerge in the academic field and, over time, are adopted by public institutions, private companies and statistical organizations, which demand their incorporation into the main data analysis programs.
This chapter presents the main procedures for the statistical analysis of continuous variables using R, appropriately incorporating the characteristics of the sampling design, such as expansion weights, stratification and conglomeration. Throughout the chapter, methods are illustrated to estimate population parameters, measure their precision, and make valid inferences from data from complex surveys. Likewise, functions and tools of the survey (Lumley, 2024) package are introduced, one of the most used resources in R for survey analysis.
Databases used in survey analysis can be found in a wide variety of formats (.xlsx, .dat, .csv, .parquet, .sas7bdat, .sav, .txt, among others). However, a good practice is to import the information from any of these formats and immediately save it in a file with the extension .rds, which is the native format of R.
.rds files allow you to store any object or information in R, such as data frames, vectors, matrices or lists, and are characterized by their complete compatibility with the R work environment. In addition, they facilitate efficient reading and writing of data, reducing loading times and avoiding compatibility problems between sessions. Another R format is .Rdata, which allows multiple objects to be saved in a single file. In contrast, .rds files store only one object, making their use more controlled and reproducible. For this reason, it is recommended to preferably work with the .rds format when documenting projects or reproducing analyses.
To illustrate the computational syntax that will be used throughout the chapter, we use the same database as in the previous chapter, which contains a sample of 2,605 records from a complex sample design. The following shows how to load a file with the .rds extension:
According to Naciones Unidas (2009, Section 7.8), it is essential that the structure of the sampling design be taken into account in the inference process when estimating official statistics based on household surveys. As emphasized in this document, ignoring this aspect can generate biased estimates and underestimated sampling errors. In this sense, statistical programs incorporate specific functionalities to handle data from surveys with complex designs (Heeringa et al., 2017, appendix A).
Once the sample of households has been loaded into R, the next step is to define the sample design from which said sample comes. For this, the srvyr (Freedman Ellis & Schneider, 2024) package will be used, which appears as a complement to the survey (Lumley, 2024) package. These libraries allow defining survey.design type objects, to which the estimation and survey analysis functions are applied, and which can be combined with the programming of the tidyverse (Wickham, Averick, et al., 2026) package. Here is an example of defining an object of type survey_design for the survey:
options(survey.lonely.UPM = "adjust")
library(survey)
library(srvyr)
survey_design <- survey_data %>%
as_survey_design(
strata = Stratum,
ids = PSU,
weights = wk,
nest = TRUE
)In this code, survey_data corresponds to the database that contains the selected sample. From this data, the function as_survey_design() converts the information into an object of class survey_design, which stores all the elements necessary to describe the sample design of the survey. These include the stratification variable (strata), which identifies the strata defined in the design, and the variable specified in ids, which represents the primary sampling units (PSUs) or clusters selected in the first stage of sampling. Likewise, the weights argument defines the variable that contains the expansion factors associated with each observation. Finally, the nest = TRUE option indicates that PSUs are nested within strata, a necessary characteristic to correctly represent the hierarchical structure of the sampling design.
## Estimation of descriptive parameters {#estimacion-de-parametros-descriptivos}
In household surveys it is common to estimate descriptive parameters associated with numerical variables, such as totals, means, ratios and standard deviations. These measures allow us to characterize different social and economic phenomena of the population, for example income, expenditure, hours worked or household size. This section presents the logic behind the syntax used for this type of estimation, along with its corresponding precision measures. ### Total estimate {#estimacion-del-total}
In R, the survey_total() function is used to estimate totals from the design object. This feature automatically incorporates expansion factors and design structure, ensuring that estimates correctly reflect the inference model.
The main arguments of the survey_total() function are x, na.rm, vartype, level, and deff. The x argument corresponds to the variable on which you want to calculate the total; It can also be left empty when seeking to obtain only the weighted total of cases in the sample. Meanwhile, na.rm indicates whether missing values (NA) should be removed before performing the calculation, with FALSE being the default value.
The vartype argument allows you to specify the type of uncertainty measure that you want to estimate, and one or more results can be requested simultaneously, such as the standard error ("se"), the confidence interval ("ci"), the variance ("var") or the coefficient of variation ("cv"). Likewise, level defines the confidence level used to calculate the confidence intervals, whose default value is 0.95. On the other hand, the confidence interval can be obtained directly by including the argument vartype = "ci" within survey_total(), or by using the function confint() on the resulting object, specifying the required confidence level. Finally, deff is a logical value that indicates whether to return the design effect (DEFF).
The following example presents how to estimate totals and their confidence intervals for different variables of interest in R, using the function survey_total(), and presented in Table 3.1:
| total | total_se | total_low | total_upp | total_deff |
|---|---|---|---|---|
| 85793667 | 4778674 | 76331414 | 95255920 | 11 |
From the table, it is obtained that the estimate of the total population for the total income is 85,793,667. This value corresponds to the point estimator obtained from the survey, incorporating the expansion factors and the characteristics of the sample design. The standard error associated with this estimate is 4,778,674, which reflects the sample variability of the estimator. From this value, a 95% confidence interval is constructed, whose lower and upper limits are 76,331,414 and 95,255,920, respectively. The design effect (deff) is equal to 11, suggesting a significant loss of accuracy mainly associated with complex design characteristics. Consequently, although the nominal sample size may be large, the effective sample size would be considerably smaller due to the correlation between observations within the clusters.
To continue illustrating the use of the survey_total() function within the srvyr approach, let us estimate total household expenditures, but now calculating the 90% confidence interval. The following code performs this estimation, presented in Table 3.2:
survey_design %>%
summarise(total = survey_total(Expenditure,
vartype = "ci", level = 0.9,
deff = TRUE))| total | total_low | total_upp | total_deff |
|---|---|---|---|
| 55677504 | 51360469 | 59994539 | 10.2 |
If the objective is to estimate total household income, but disaggregated by sex, the function group_by() can be used to group by the variable of interest, together with cascade() from the srvyr library, which allows adding a row with the general total at the end of the table. In this way, category totals and the overall total can be easily obtained within the same summary framework, as shown in Table 3.3.
survey_design %>%
group_by(Sex) %>%
cascade(total = survey_total(Income,
vartype = c("se", "ci"),
level = 0.95),
.fill = "Total income")| Sex | total | total_se | total_low | total_upp |
|---|---|---|---|---|
| Female | 44153820 | 2324452 | 39551172 | 48756467 |
| Male | 41639847 | 2870194 | 35956576 | 47323118 |
| Total income | 85793667 | 4778674 | 76331414 | 95255920 |
As seen in the previous codes, a practical way to obtain the estimates of the total, its standard error and its confidence interval is through the argument vartype, specifying the options "se" and "ci" respectively.
### Estimation of the mean {#estimacion-de-la-media}
The estimation of the population average occupies a central place in household surveys, since many of the indicators used for socioeconomic analysis correspond to average values, such as average household income, per capita spending or average hours worked. This type of parameters allows us to summarize the general behavior of the variables of interest and describe the central tendencies of the target population. According to Gutiérrez (2016), the population mean estimator can be expressed as a non-linear ratio between two population totals, which are estimated as follows: \[ \hat{\bar{y}}= \frac{\hat{t}_y}{\hat{N}} = \frac{\sum_{h}\sum_{i}\sum_{k} w_{h i k} \ y_{hik}}{\sum_{h}\sum_{i}\sum_{k} w_{h i k}} \]
Since \(\hat{\bar{y}}\) is not a linear statistic, there is no closed formula for the variance of this estimator. For this reason, it is necessary to resort to resampling methods or Taylor linearization to approximate its variance. In this particular case, using Taylor series, the variance is defined as: \[ \widehat{Var}\left(\hat{\bar{y}}\right) \approx \frac{\widehat{Var}\left(\hat{t}_y\right)+\hat{\bar{y}}^{2}\ \widehat{Var}\left(\hat{N}\right)-2\ \hat{\bar{y}} \ \widehat{Cov}\left(\hat{t}_y,\hat{N}\right)}{\hat{N}^{2}} \]
As can be seen, the calculation of the variance of the population mean involves complex analytical components, such as the covariance between the estimated total and the estimated population size. However, the srvyr package in R facilitates these calculations by incorporating functions that directly estimate the mean, its standard error, the confidence interval, and the design effect. Below is the syntax to estimate the average household income, along with its 95% confidence interval, presented in Table 3.4:
survey_design %>%
summarise(mean = survey_mean(Income,
vartype = c("se", "ci"),
level = 0.95,
deff = TRUE))| mean | mean_se | mean_low | mean_upp | mean_deff |
|---|---|---|---|---|
| 571 | 28.5 | 515 | 627 | 8.82 |
As can be seen, the arguments used in the survey_mean() function are similar to those of survey_total(). In this case, vartype allows you to obtain the standard error ("se") and the confidence interval ("ci"), while the argument deff = TRUE requests the calculation of the design effect. Similarly, it is possible to estimate the average household expenses, using the same code structure, presented in Table 3.5:
survey_design %>%
summarise(mean = survey_mean(Expenditure,
vartype = c("se", "ci"),
level = 0.95,
deff = TRUE))| mean | mean_se | mean_low | mean_upp | mean_deff |
|---|---|---|---|---|
| 371 | 13.3 | 344 | 397 | 6.02 |
Estimates of the mean by subgroups can also be made following the same scheme shown previously. Particularly, household expenses discriminated by sex are presented in Table 3.6:
survey_design %>%
group_by(Sex) %>%
cascade(mean = survey_mean(Expenditure,
level = 0.95,
vartype = c("se", "ci")),
.fill = "Mean expenditure") %>%
arrange(desc(Sex))| Sex | mean | mean_se | mean_low | mean_upp |
|---|---|---|---|---|
| Mean expenditure | 371 | 13.3 | 344 | 397 |
| Male | 374 | 16.1 | 343 | 406 |
| Female | 367 | 12.3 | 343 | 391 |
In an analogous manner, the average expenditure by area can be calculated, which is presented in Table 3.7:
survey_design %>%
group_by(Zone) %>%
cascade(mean = survey_mean(Expenditure,
level = 0.95,
vartype = c("se", "ci")),
.fill = "Mean expenditure") %>%
arrange(desc(Zone))| Zone | mean | mean_se | mean_low | mean_upp |
|---|---|---|---|---|
| Urban | 460 | 22.2 | 416 | 504 |
| Rural | 274 | 10.3 | 254 | 294 |
| Mean expenditure | 371 | 13.3 | 344 | 397 |
Likewise, it is possible to make a combined estimate by sex and area, presented in Table 3.8:
survey_design %>%
group_by(Zone, Sex) %>%
cascade(mean = survey_mean(Expenditure,
level = 0.95,
vartype = c("se", "ci")),
.fill = "Mean expenditure") %>%
arrange(desc(Zone), desc(Sex))| Zone | Sex | mean | mean_se | mean_low | mean_upp |
|---|---|---|---|---|---|
| Urban | Mean expenditure | 460 | 22.2 | 416 | 504 |
| Urban | Male | 470 | 27.0 | 416 | 523 |
| Urban | Female | 451 | 20.1 | 411 | 491 |
| Rural | Mean expenditure | 274 | 10.3 | 254 | 294 |
| Rural | Male | 275 | 10.2 | 255 | 296 |
| Rural | Female | 273 | 11.6 | 250 | 296 |
| Mean expenditure | Mean expenditure | 371 | 13.3 | 344 | 397 |
3.0.1 Ratio estimation
A particular case of nonlinear functions of totals is the population ratio, defined as the quotient between two population totals associated with characteristics of interest. Ratios allow the relationship between two variables to be expressed, which is especially useful for building comparative and monitoring indicators.
This type of parameter is widely used in household surveys to construct relative indicators, for example, the number of men for each woman, the proportion of employed people compared to the working-age population or the average number of pets per household. Ratios are also used in international frameworks. For example, Indicator 2.1.1 of the Sustainable Development Goals (SDGs) (prevalence of undernourishment) is calculated from the ratio between food consumption, measured in calories ingested, and the minimum energy requirements of the diet, determined according to age, sex and level of physical activity.
Since the ratio corresponds to the quotient between two unknown population parameters, both the numerator and the denominator must be estimated from the sample (Bautista, 1998). Being \(Y\) the total of the variable \(y\) and \(X\) the total of the variable \(x\), the population ratio is defined as \(R = \frac{t_y}{t_x}\) and its point estimator in the framework of complex sampling is expressed as: \[ \hat{R} = \frac{\hat{t}_y}{\hat{t}_x} = \frac{\sum_{h}\sum_{i}\sum_{k} w_{hik} \ y_{hik}} {\sum_{h}\sum_{i}\sum_{k} w_{hik} \ x_{hik}} \]
However, since \(\hat{R}\) is a quotient between two estimators, that is, two random variables, the calculation of its variance is not trivial. For this, Taylor linearization is used, as shown by Gutiérrez (2016), or resampling methods. In particular, the survey_ratio function implements the estimation of ratios and their variances in the framework of complex surveys.
A direct example in household surveys is the expenditure/income ratio, which helps identify consumption patterns and levels of economic sustainability of households. The following illustrates how to estimate the ratio between household expenditure and household income, presented in Table 3.9:
survey_design %>%
summarise(ratio = survey_ratio(numerator = Expenditure,
denominator = Income,
level = 0.95,
vartype = c("se", "ci"))
)| ratio | ratio_se | ratio_low | ratio_upp |
|---|---|---|---|
| 0.649 | 0.023 | 0.603 | 0.695 |
In this case, the expenditure variable was specified as the numerator (numerator), the income variable as the denominator (denominator), the confidence level (level) used to construct the confidence intervals, and the precision measures required by the argument (vartype). As can be seen, the ratio between expenditure and income is, approximately, 0.65. Which implies that for every 100 monetary units that enter the home, 65 units are spent, achieving a 95% confidence interval of 0.60 and 0.69.
If now the objective is to estimate the ratio between women and men in the example base, it is done as follows, presented in Table 3.10:
survey_design %>%
summarise(ratio = survey_ratio(numerator = (Sex == "Female"),
denominator = (Sex == "Male"),
level = 0.95,
vartype = c("se", "ci"))
)| ratio | ratio_se | ratio_low | ratio_upp |
|---|---|---|---|
| 1.11 | 0.035 | 1.04 | 1.18 |
Given that the sex variable in the database is categorical, it was necessary to use indicator variables to calculate the ratio, using Sex == "Female" to identify women and Sex == "Male" for men. The results show that in the analyzed population there are a greater number of women than men, obtaining an estimated ratio of 1.11. This indicates that, for every 100 men, there are approximately 111 women. Furthermore, the 95% confidence interval suggests that this ratio could vary between 1.05 and 1.18. If you want to make the ratio of women and men but in rural areas, it would be done as follows, presented in Table 3.11:
rural_subset <- survey_design %>%
filter(Zone == "Rural")
rural_subset %>%
summarise(ratio = survey_ratio(numerator = (Sex == "Female"),
denominator = (Sex == "Male"),
level = 0.95,
vartype = c("se", "ci"))
)| ratio | ratio_se | ratio_low | ratio_upp |
|---|---|---|---|
| 1.07 | 0.035 | 0.997 | 1.14 |
Now, another analysis of interest is to estimate the expense ratio but only in the female population. The computational codes are presented below, as shown in Table 3.12.
female_subset <- survey_design %>%
filter(Sex == "Female")
female_subset %>%
summarise(ratio = survey_ratio(numerator = Expenditure,
denominator = Income,
level = 0.95,
vartype = c("se", "ci"))
)| ratio | ratio_se | ratio_low | ratio_upp |
|---|---|---|---|
| 0.658 | 0.02 | 0.619 | 0.698 |
The result is that for every 100 monetary units that women receive, 66 are spent with a confidence interval between 0.62 and 0.69. Finally, similarly for men, the expense ratio is very similar to that for women, as shown in Table 3.13.
male_subset <- survey_design %>%
filter(Sex == "Male")
male_subset %>%
summarise(ratio = survey_ratio(numerator = Expenditure,
denominator = Income,
level = 0.95,
vartype = c("se", "ci"))
)| ratio | ratio_se | ratio_low | ratio_upp |
|---|---|---|---|
| 0.639 | 0.029 | 0.582 | 0.696 |