The R Package smicd: Statistical Methods for Interval-Censored Data

The package allows the use of two new statistical methods for the analysis of intervalcensored data: 1) direct estimation/prediction of statistical indicators, and 2) linear (mixed) regression analysis. Direct estimation of statistical indicators, for instance poverty and inequality indicators, is facilitated by a non-parametric kernel density algorithm. The algorithm is able to account for weights in the estimation of statistical indicators. The standard errors of the statistical indicators are estimated with a non-parametric bootstrap. Furthermore, the package offers statistical methods for the estimation of linear and linear mixed regression models with an interval-censored dependent variable, particularly random slope and random intercept models. Parameter estimates are obtained through a stochastic expectation-maximization algorithm. Standard errors are estimated using a non-parametric bootstrap in the linear regression model and by a parametric bootstrap in the linear mixed regression model. To handle departures from the model assumptions, fixed (logarithmic) and data-driven (Box-Cox) transformations are incorporated into the algorithm.


Introduction
Interval-censored or grouped data occurs when only the lower A k−1 and upper A k interval bounds (A k−1 , A k ) of a variable are observed and its true value remains unknown. Instead of measuring the variable of interest on a continuous scale, for instance income data, the scale is divided into n k intervals. The variable k (1 ≤ k ≤ n k ) indicates in which of the n k intervals an observation falls into. This leads to a loss of information since the shape of the distribution within the intervals remains unknown. In the field of survey statistics, asking for interval-censored data is often done in order to avoid item non-response and thus increase data quality. Item non-response is avoided because interval-censored data offers a higher level of data privacy protection (Hagenaars and Vos, 1988;Moore and Welniak, 2000). Among others, popular surveys and censuses that collect intervalcensored data are the German Microcensus (Statistisches Bundesamt, 2017), the Colombian census (Departamento Administrativo Nacional De Estadística, 2005) and the Australian census (Australian Bureau of Statistics, 2011). While item non-response is reduced or avoided, the statistical analysis of the data requires more elaborate mathematical methods. Even statistical indicators that are easily calculated for metric data, e.g., the mean, cannot be estimated using standard formulas (Fahrmeir et al., 2016). Also, estimating linear and linear mixed regression models which are applied in many fields of statistics requires advanced statistical methods when the dependent variable is interval censored. Therefore, the presented R package implements three major functions: kdeAlgo() to estimate statistical indicators (e.g., the mean) from interval-censored data, semLm() and semLme() to estimate linear and linear mixed regression models with an interval-censored dependent variable. The package code and the open source contribution guidelines for the package are available on GitHub. Potential code contributions, feature requests and bugs can be reported there by creating issues.
For the estimation of statistical indicators from interval-censored data different approaches are described in the literature. These approaches can be broadly categorized into four groups: Estimation on the midpoints (Fahrmeir et al., 2016), linear interpolation of the distribution function, non-parametric modeling via splines (Berger and Escobar, 2016) and fitting a parametric distribution function to the censored data (Dagum, 2008;McDonald, 1984;Bandourian et al., 2002). Some of these methods are implemented in R packages available on the Comprehensive R Archive Network (CRAN). The method of linear interpolation is implemented for the estimation of quantiles in the R package actuar (Goulet et al., 2020;Dutang et al., 2008). The package also enables the estimation of the mean on the interval midpoints. Fitting a parametric distribution to interval-censored data can be done by using the R package fitdistrplus (Delignette-Muller et al., 2020;Delignette-Muller and Dutang, 2015).
In survey statistics, interval-censored data is often collected for income or wealth variables. Thus, the performance of the above-mentioned methods is commonly evaluated by simulation studies that rely on data that follows some kind of income distribution. The German statistical office (DESTATIS) uses the method of linear interpolation for the estimation of statistical indicators from interval-censored income data collected by the German Microcensus. This approach gives the same results as assuming a uniform distribution within the income intervals. Estimation results are reasonably accurate if the estimated indicators do not depend on the whole shape of the distribution, e.g., the median (Lenau and Münnich, 2016). Fitting a parametric distribution to the data enables the estimation of indicators that rely on the whole shape of the distribution. This method works well when the data is censored to only a few equidistant intervals (Lenau and Münnich, 2016). Non-parametric modeling via splines shows especially good results for a high number of intervals in ascending order (Lenau and Münnich, 2016). However, according to Lenau and Münnich (2016) all of the above-mentioned methods show large biases and variances when the estimation is based on a small number of intervals. Therefore, a novel kernel density estimation (KDE) algorithm is implemented in the smicd package that overcomes the drawbacks of the previously mentioned methods (Walter, 2019(Walter, , 2020. The algorithm bases the estimation of statistical indicators on pseudo samples that are drawn from a fitted non-parametric distribution. The method automatically adapts to the shape of the true unknown distribution and provides reliable estimates for different interval-censoring scenarios. It can be applied via the function kdeAlgo().
Similarly to the direct estimation of statistical indicators from interval-censored data, there exists a variety of ad-hoc approaches and explicitly formulated mathematical methods for the estimation of linear regression models with an interval-censored dependent variable. The following methods and approaches are used for handling interval-censored dependent variables within linear regression models: Ordinary least squares (OLS) regression on the midpoints (Thompson andNelson, 2003), ordered logit-or probit-regression (McCullagh, 1980), and regression methodology formulated for left-, right-, and interval-censored data (Tobin, 1958;Rosett and Nelson, 1975;Stewart, 1983). All of these methods are implemented in different R packages available on CRAN. OLS regression on the midpoints is applicable by using the lm() function from the stats Package (R Core Team, 2020), ordered logit regression is implemented in the MASS package (Ripley, 2019;Venables and Ripley, 2002), and interval regression is implemented in the survival (Therneau, 2020;Therneau and Grambsch, 2000) package.
While OLS regression on the midpoints of the intervals is easily applied, it comes with the disadvantage of giving biased estimation results (Cameron, 1987). This approach disregards the uncertainty stemming from the unknown true distribution of the data within the intervals and therefore leads to biased parameter estimates. Its performance relies on the number of intervals and estimation results are only comparable to more advanced methods when the number of intervals is very large (Fryer and Pethybridge, 1972). Conceptualizing the model as an ordered logit or probit regression is feasible by treating the dependent variable as an ordered factor variable (McCullagh, 1980). However, this approach also neglects the unknown distribution of the data within the intervals. Furthermore, the predicted values are not on a continuous scale but are in terms of probability of belonging to a certain group. To overcome these disadvantages and obtain unbiased estimation results Stewart (1983) introduces regression methodology for models with an interval-censored dependent variable. Walter (2019) further develop his approach and introduce a novel stochastic expectationmaximization (SEM) algorithm for the estimation of linear regression models with an interval-censored dependent variable that is implemented in the smicd package. The model parameters are unbiasedly estimated as long as the model assumptions are fulfilled. The function semLm() provides the SEM algorithm and enables the use of fixed (logarithmic) and data-driven (Box-Cox) transformations (Box and Cox, 1964). The Box-Cox transformation automatically adapts to the shape of the data and transforms the dependent variable in order to meet the model assumption.
In order to analyze longitudinal or clustered data (e.g., students within schools) linear mixed regression models are applicable. These kinds, of models control for the correlated structure of the data by including random effects in addition to the usual fixed effects. In order to deal with an interval-censored dependent variable in linear mixed regression models there are several approaches described in the literature. Linear mixed regression models, just like linear regression models, can be estimated on the interval midpoints of the censored-dependent variable. Furthermore, conceptualizing the model as an ordered logit or probit regression model is feasible (Agresti, 2010). These approaches inherit the same advantages and disadvantages as previously discussed. Linear mixed regression on the midpoints can be applied by the lme4 (Bates et al., 2020b(Bates et al., , 2015 or nlme (Pinheiro et al., 2020) package and the ordered logit regression is implemented in the ordinal package (Christensen, 2019). To my knowledge, there are no R packages for the estimation of linear mixed regression models with an interval-censored dependent variable. Therefore, the package smicd contains the SEM algorithm proposed by Walter (2019) for the estimation of linear mixed regression models with an intervalcensored dependent variable. If the model assumptions are fulfilled, the method gives unbiased estimation results. The function semLme() enables the estimation of the regression parameters and it also allows for the usage of the logarithmic and Box-Cox transformation in order to fulfill the model assumptions (Gurka et al., 2006).
The paper is structured into two main sections. Section 2.2 deals with the direct estimation of statistical indicators from interval-censored data whereas Section 2.3 introduces linear and linear mixed regression models with an interval-censored dependent variable. Both sections have been divided into three subsections: first the statistical methodology is introduced, then the core functions of the smicd package are presented, and finally, illustrative examples with two different data sets are provided. In Section 2.4 the main results are summarized and an outlook is given.

Direct estimation of statistical indicators
In the following three subsections, the methodology for the direct estimation of statistical indicators from interval-censored data is introduced, the core functionality of the function kdeAlgo() is presented and statistical indicators are estimated using the synthetic EU-SILC (European Union Statistics on Income and Living Conditions) data set from Austria.

Methodology: Direct estimation of statistical indicators
In order to estimate statistical indicators from interval-censored data the proposed algorithm generates metric pseudo samples of an interval-censored variable. These pseudo samples can be used to estimate any statistical indicator. They are drawn from a non-parametrically estimated kernel density. Kernel density estimation was first introduced by (Rosenblatt, 1956) and (Parzen, 1962). By its application the density f (x) of a continuous independently and identically distributed random variable is estimated without assuming any distributional shape of the data. The estimator is defined aŝ where K (·) is a kernel function, h > 0 the bandwidth and x = {x 1 , x 2 , . . . , x n } denotes a sample of size n. The performance of the estimator is determined by the optimal choice of h. The selection of an optimal h is widely discussed in the literature, see Jones et al. (1996);Loader (1999); Zambom and Dias (2012). When working with interval-censored data, a standard KDE cannot be applied since x is not observed on a continuous scale. Nevertheless, its unobserved true distribution is of continuous form. As an ad hoc solution the densityf h (x) can be estimated based on the interval midpoints. The resulting density estimate will be spiky unless the bandwidth is sufficiently large. A large bandwidth, however, leads to a loss of information (Wang and Wertelecki, 2013). Therefore, Walter (2019) proposes an iterative KDE algorithm for density estimation from interval-censored data. The approach is based on Groß et al. (2017) who introduce a similar KDE algorithm in a two-dimensional setting with an equidistant interval width. Walter (2019) shows that the algorithm can be adjusted to one-dimensional data with an arbitrary class width. For the estimation of linear and non-linear statistical indicators the unknown distribution of x has to be reconstructed by using the observed interval k = {k 1 , k 2 , . . . , k n } that an observation falls into. From Bayes' theorem (Bayes, 1763) it follows that the conditional distribution of (x|k) is: Since π (x) is unknown it is replaced by a kernel density estimatef h (x).

Estimation and computational details
To fit the model, pseudosamples of x i are drawn from the conditional distribution where I (·) denotes the indicator function. The conditional distribution of π (x i |k i ) is given by the product of a uniform distribution and density f (x i ). As the density is unknown it is replaced by an estimatef h (x), which is obtained by the KDE. In particular, x i is repeatedly drawn from the given interval A k i −1 , A k i by using the current density estimatef h (x) as a sampling weight. The explicit steps of the iterative algorithm as given in Walter (2019) are stated below: 1. Use the midpoints of the intervals as pseudox i for the unknown x i . Estimate a pilot estimate off h (x), by applying KDE. Note: Choose a sufficiently large bandwidth h in order to avoid rounding spikes.
2. Evaluatef h (x) on an equal-spaced grid G = {g 1 , . . . , g j } with grid points g 1 , . . . , g j . The width of the grid is denoted by δ g . It is given by and the grid is defined as: For open-ended intervals, e.g., (15000, +∞) the upper bound has to be replaced by a finite number.
Walter (2019) shows through model-based simulations that a value of three times the value of the lower bound (15000, 45000) gives appropriate estimation results when working with income data. The variance of the statistical indicators is estimated by bootstrapping. Bootstrap methods were first introduced by Efron (1979). These methods serve as an estimation procedure when the variance cannot be stated as closed-form solution (Shao and Tu, 1995). While bootstrapping avoids the problem of the non-availability of a closed-form solution, it comes with the disadvantage of long computational times. In the package, a non-parametric bootstrap that accounts for the additional uncertainty coming from the interval-censored data is implemented. This non-parametric bootstrap is introduced in Walter (2019).

Core functionality: Direct estimation of statistical indicators
The presented KDE algorithm is implemented in the function kdeAlgo() (see Table 1). The arguments and default settings of kdeAlgo() are briefly summarized in Table 2. The function gives back an S3 object of class "kdeAlgo". A detailed explanation of all components of an "kdeAlgo" object can be found in the package documentation. The generic functions plot() and print() can be applied to "kdeAlgo" objects to output the main estimation results (see Table 1). In the next section the function kdeAlgo() is used to estimate a variety of statistical indicators from interval-censored EU-SILC data and its arguments are explained in more detail.

kdeAlgo()
Estimates statistical indicators and its standard errors from interval-censored data plot() Plots convergence of the estimated statistical indicators and estimated density of the pseudox i print() Prints estimated statistical indicators and its standard errors

Example: Direct estimation of statistical indicators
To demonstrate the function kdeAlgo(), the equivalized household income and the corresponding household sample weight from the Austrian synthetic EU-SILC survey data set is used. The data set is included in the laeken package (Alfons et al., 2020;Alfons and Templ, 2013). Its synthetic nature makes it unusable for inferential statistics. However, the data set has the advantage over the scientific use file by being freely available which enables the easy reproducibility of the presented example. Since the total disposable household income is measured on a continuous scale, it is censored to 22 intervals for demonstration purposes. For a realistic censoring scheme the interval bounds are chosen such that they closely follow the interval bounds used in the German Microcensus from 2013 (Statistisches Bundesamt, 2014). The German Microcensus is a representative household survey that covers 830,000 persons in 370,000 households (1% of the German population) in which income is only collected as interval-censored variable (Statistisches Bundesamt, 2016).
In a first step the variable total disposable household income called hhincome_net is interval censored according to 22 intervals using the function cut(). The vector of interval bounds is called intervals and the newly obtained interval-censored income variable is called c.hhincome. The variable c.hhincome is assigned to the argument xclass and the vector of interval bounds intervals is assigned to the argument classes. The default settings of the arguments burnin, samples, bw, evalpoints, adjust and upper are retained. Simulation results from Walter (2019) and Groß et al. (2017) show that these settings give good results when working with income data. Changing these arguments has an impact on the performance of the KDE algorithm. As default, the statistical indicators: mean, Gini coefficient, headcount ratio (HCR), the quantiles (10%, 25%, 50%, 75%, 90%), the poverty gap (PGAP) and the quintile share ratio (QSR) are estimated (Gini, 1912;Foster et al., 1984). The HCR and PGAP rely on a poverty threshold. The default choice of the threshold argument is 60% of the median of the target variable as suggested by Eurostat (2014). Besides the mentioned indicators, any other statistical indicator can be estimated via the argument custom_indicator. In the example the argument is assigned a list that holds functions to estimate the 5% and 95% quantile. The custom indicators must depend on the target variable, the threshold (even if it is not needed for the specified indicator) and optionally on the weights argument, if the estimation of a weighted indicator is required. To estimate the standard errors of all indicators bootstrap.se = TRUE and the number of bootstrap samples is 100 (the default value as suggested in Walter (2019)). Lastly, the household weight (hhweight) is assigned to the argument weights in order to estimate weighted statistical indicators. It can also be controlled for households of different sizes, by assigning oecd a variable with household equivalence weights. By applying the print() function to the "kdeAlgo" object the estimated statistical indicators (default and custom indicators) as well as their standard errors are printed. For instance, in this example the estimated mean is about 1,658 Euro and its standard error is 8.486. For demonstration purposes, the statistical indicators are also estimated using the continuous household income variable from the synthetic EU-SILC data set (Table 3). The estimation results of the KDE algorithm using the interval-censored data are very close to those based on the continuous data. Slightly larger deviations are observable for the more extreme quantiles. This is due to the fact that these quantile estimates fall into intervals with a lower number of observations (compared to the other quantile estimates). Estimation results for these quantiles could potentially be further improved by increasing the number of evalpoints of the kdeAlgo(). In Walter (2019) the performance of the KDE algorithm is evaluated via detailed simulation studies. By applying the function plot(), "kdeAlgo" objects can be plotted. Thereby, convergence plots for all estimated statistical indicators and a plot of the estimated final density are obtained.

R> plot(Indicators)
The R Journal Vol. XX/YY, AAAA 20ZZ ISSN 2073-4859 Figure 1 shows convergence plots for two of the estimated indicators. Additionally, a plot of the estimated final density with a histogram of the observed data in the background is given in Figure 2. In Figure 1 the estimated statistical indicator (Gini, 10% quantile) for each iteration step of the KDE algorithm and the average over the estimates up to iteration step M (excluding the burn-in iterations) are plotted. A vertical line marks the end of the burn-in period. The horizontal line gives the value of the final estimate (average over the M iterations). All convergence plots indicate that the number of iterations is chosen sufficiently large for the estimates to converge.
If convergence were not achieved the arguments burnin and samples should be increased. It is notable that the estimated 10% quantile has the same value for almost all iterations steps. This is the case because the quantile, as any other statistical indicator, is estimated using the pseudo samples that are drawn on 4,000 grid points G. Estimating a quantile on only 4,000 unique outcomes (pseudo values) leads to equal quantile estimates for numerous iteration steps of the KDE algorithm. Lastly, it should be mentioned that the computation time can be very long if the estimation of the standard errors is enabled. Hence, if the estimation of the standard errors is not required, the argument bootstrap.se should be set to FALSE. Furthermore, it should always be checked how many iterations are needed for the desired statistical indicator to converge. Reducing the number of required iterations (arguments burnin and samples) lowers the computation time significantly (with and without standard errors).

Regression analysis
In the following three subsections the statistical methodology for linear and linear mixed regression models with an interval-censored dependent variable is introduced, the core functionality of the functions semLM() and semLME() is presented and examination scores of students from schools in London are exemplary modeled.

Methodology: Regression analysis
The theoretical introduction of the new regression method, proposed by Walter (2019), is presented for linear mixed regression models. The theory for linear regression models can be obtained by simplifying the introduced method. In its standard form the linear mixed regression model serves to analyze the linear relationship between a continuous dependent variable and some independent variables (Goldstein, 2010). Random parameters (random slopes and random intercepts) are included in the model to account for correlated data e.g., students within schools. The model in matrix notation (Laird and Ware, 1983) is given by where y is a n × 1 column vector of the dependent variable, n is the sample size, X is a n × p matrix where p is equal to the number of predictors, β is a column vector of the fixed effects regression parameters of size p × 1, Z is the n × q design matrix with q random effects, v is a q × 1 vector of random effects, and e is the residual vector of size n × 1. The distribution of the random effects is given by , and the distribution of the residuals is given by e ∼ N (0, R) with R = I n σ 2 e where I n is the identity matrix and σ 2 e is the residual variance. The random effects v and the residuals e are assumed to be independent. For a more detailed introduction of mixed models see Searle et al. (1992);McCulloch et al. (2008); Snijders and Bosker (2011). In the case of an interval-censored dependent variable the parameters of Model (1) have to be estimated without observing y on a continuous scale. Instead, only the interval identifier k, now defined as n × 1 column vector, is observed. Open-ended interval bounds A 0 = −∞ and A n k = +∞ and unequal interval widths are allowed. Since the true distribution of y is unknown the aim is to reconstruct the distribution of y using the known intervals k and the linear relationship stated in Model (1). As presented in Walter (2019) in order to reconstruct the unknown distribution of f (y|X, Z, v, k, θ), where θ = (β, R, G), the Bayes theorem (Bayes, 1763) is applied.
with f (k|y, X, Z, v, θ) = f (k|y) because the conditional distribution of the interval identifier k only depends on y. It is given by f (k|y) = r with r being a n × 1 column vector r = (r 1 , r 2 , . . . , r n ) T with The relationship in Equation (2) follows from the linear mixed model assumptions (Model (1)). The unknown parameters θ = (β, R, G) are estimated based on pseudo samplesỹ (since y is unknown) that are iteratively drawn from f (y|X, Z, v, k, θ). The next subsection states the computational details of the SEM algorithm.

Estimation and computational details
To fit Model (1), the parameter vectorθ = β, R, G is estimated and pseudo samples of the unknown y are iteratively generated by the following SEM algorithm. The pseudo samplesỹ are drawn from the conditional distribution where I (·) denotes the indicator function. Hence, for y with explanatory variables X the corresponding y is drawn from N (Xβ + Zv, R) conditional on the given interval (A k−1 ≤ y ≤ A k ). Ifθ is estimated the conditional distribution f (y|X, Z, v, k, θ) follows a two-sided truncated normal distribution. Its probability density function equalŝ withμ = X β + Zv. φ (·) denotes the probability density function of the standard normal distribution and Φ(·) denotes its cumulative distribution function. From its definition it follows that Φ The steps of the SEM algorithm as described in Walter (2019) are: 1. Estimateθ = β, R, G from Model (1) using the midpoints of the intervals as substitutes for the unknown y. The parameters are estimated by restricted maximum likelihood theory (REML) (Thompson, 1962).

Stochastic step:
For i = 1, . . . , n, draw randomly from N X β + Zv, R within the given interval (A k−1 ≤ y ≤ A k ) (the two-sided truncated normal distribution given in Equation (3)) obtaining (ỹ, X, Z). The drawn pseudoỹ are used as replacements for the unobserved y.

Maximization step:
Re-estimate the parameter vectorθ from Model (1) by using the pseudo samples (ỹ, X, Z) from Step 2. Again, parameter estimation is carried out by REML. If open-ended intervals A 0 = −∞ and A n k = +∞ are present, the midpoints M 1 and M n k of these intervals in iteration

Iterate
Step 1 are computed as follows: The R Journal Vol. XX/YY, AAAA 20ZZ ISSN 2073-4859 These midpoints serve as proxies for the unknown interval midpoints in Step 1 of the algorithm. The SEM algorithm for the linear regression model is obtained by simplifying the conditional distribution f (y|X, Z, v, θ) ∼ N (Xβ + Zv, R) to f y|X, β, σ 2 e ∼ N Xβ, σ 2 e according to the model assumptions of a linear regression model. In the SEM algorithm for linear models it is then drawn from N Xβ, σ 2 e within the given interval.
The standard errors of the regression parameters are estimated using bootstrap methods. For the linear regression model a non-parametric bootstrap (Efron and Stein, 1981;Efron, 1982;Tibshirani, 1986, 1993) and for the linear mixed regression model a parametric bootstrap (Wang et al., 2006;Thai et al., 2013) is used to estimate the standard errors. The non-parametric as well as the parametric bootstrap are further developed to account for the additional uncertainty that is due to the interval-censored dependent variable. Both newly proposed bootstraps are available in the smicd package and the theory is explained in (Walter, 2019).
To assure that the model assumptions are fulfilled the logarithmic and the Box-Cox transformations are incorporated into the function semLm() and semLme().

Core functionality: Regression analysis
The introduced SEM algorithm is implemented in the functions described in Table 4. The arguments and default settings of the estimation functions semLm() and semLme() are summarized in Table 5. Both functions return a an S3 object of class "sem" "lm" or "sem" "lme". A detailed explanation of all the components of these objects can be found in the smicd package documentation. The generic functions plot(), print() and summary() can be applied to objects of class "sem" "lm" and "sem" "lme" in order to summarize the main estimation results. In the next section the functionality of semLm() and semLme() is demonstrated based on an illustrative example.

semLm()
Estimates linear regression models with an interval-censored dependent variable semLme() Estimates linear mixed regression models with an interval-censored dependent variable plot() Plots convergence of the estimated parameters and estimated density of the pseudoỹ from the last iteration step print() Prints basic information of the estimated linear and linear mixed regression models summary() Summary of the estimated linear and linear mixed regression models

Example: Regression analysis
To demonstrate the functions semLm() and semLme() the famous London school data set that is analyzed in Goldstein et al. (1993) is used. The data set contains the examination results of 4,059 students from 65 schools in six Inner London Education Authorities. The data set is available in the R package mlmRev (Bates et al., 2020a) and also included in the package smicd. The variables used in the following example are: general certificate of secondary examination scores (examsc), the standardized London reading test scores at the age of 11 years (standLRT), the sex of the student (sex), and the school identifier (school). In the original data set the variable examsc is measured on a continuous scale. In order to demonstrate the functionality of the functions semLm() and semLme() the variable is arbitrarily censored to nine intervals. As before, the censoring is carried out by the function cut() and the vector of interval bounds is called intervals.
R> intervals 1.5,2.5,3.5,4.5,5.5,6.5,7.7,8.5,Inf) R> Exam$examsc.class <-cut(Exam$examsc, intervals) The newly created interval-censored variable is called examsc.class. The distribution is visualized by applying the function The formula argument is assigned the model equation, where examsc.class is regressed on standLRT and sex. The argument data is assigned the name of the data set Exam and the vector of interval bounds intervals is assigned to the classes argument. The arguments burnin and samples are left as defaults. The specified number of default iterations is sufficiently large for most regression models, however, convergence of the parameters has to be checked by plotting the estimation results with the function plot() after the estimation. No transformation is specified for the interval-censored dependent variable therefore trafo is assigned its default value. The argument adjust is only relevant if the Box-Cox transformation trafo="bc" is chosen. In this case the number of iterations for the estimation of the Box-Cox transformation parameter λ can be specified by this argument. The convergence of the transformation parameter λ has to be checked using the function plot(). More information on the Box-Cox transformation and on the estimation of the transformation parameter is given in Walter (2019). For the estimation of the standard errors of the regression parameters the argument bootstrap.se is set to TRUE. The number of bootstrap samples b is 100, its default value, which again is reasonable for most settings. A summary of the estimation results is obtained by the application of the function summary( The output shows the function call, the estimated regression coefficients, the bootstrapped standard errors, and the confidence intervals as well as the R-squared and the adjusted R-squared. Furthermore, the output reminds the user that the dependent variable is censored to nine intervals. All estimates are interpreted as in a linear regression model with a continuous dependent variable, hence, if standLRT increases by one unit and all other parameters are kept constant, examsc.class increases by 0.59 on average. The bootstraped confidence intervals indicate that all regressors have a significant effect on the dependent variable. By using the generic function plot() on an object of class "sem" "lm" convergence plots of each estimated regression parameter and of the estimated residual variance are obtained. Furthermore, the density of the generated pseudoỹ variable from the last iteration step is plotted with a histogram of the observed distribution of the interval-censored variable examsc.class in the background.  In Figure 3 a selection of convergence plots is given and in Figure 4 the density of the pseudoỹ from the last iteration step of the SEM algorithm is plotted. In the convergence plots the estimated parameter and the average up to iteration step M (excluding B) are plotted for each iteration step of the SEM algorithm. A vertical line indicates the end of the burn-in period (40 iterations). The final parameter estimate marked by the horizontal line is obtained by averaging the M (SEM) additional iterations (200). The selected 240 iterations are enough to obtain reliable estimates in this example, because the estimates have converged.

R> plot(LM)
As already mentioned, the smicd package also enables the estimation of linear mixed regression models by the function semLme(). In the London school data set students are nested within schools, therefore it is necessary to control for the correlation within-schools. In order to do that the variable school is specified as a random intercept. If necessary, a random slope parameter could also be included in the model. Again, the variable sex is included as an additional regressor. Hence, the formula argument is assigned the following model equation examsc.class ∼ standLRT + sex + (1|school). So far, the function semLme() enables the estimation of linear mixed models with a maximum of one random slope and one random intercept parameter. Regarding all other arguments, the same specifications as before are made.
R> LME <-semLme( + formula = examsc.class~standLRT + sex + (1|school), + data = Exam, classes = intervals, bootstrap.se = TRUE + ) By using the generic function summary() the estimation results are printed. In addition to the fixed effects, the estimated random effects are obtained as in the lme4 and nlme packages. Since the R-squared and the adjusted R-squared are not defined for mixed models the summary() function prints the marginal R-squared and conditional R-squared (Nakagawa and Schielzeth, 2013;Johnson, 2014 Again, interpretation is the same as in linear mixed models with a continuous dependent variable. By applying the generic function plot() to an "sem" "lme" object the same plots as for the linear regression model are plotted.

Discussion and outlook
Asking for interval-censored data can lead to lower item non-response rates and increased data quality. While item non-response is potentially avoided, applying traditional statistical methods becomes infeasible because the true distribution of the data within each interval is unknown. The functions of the smicd package enable researchers to easily analyze this kind of data. The paper briefly introduces the new statistical methodology and presents, in detail, the core functions of the package: • kdeAlgo() for the direct estimation of any statistical indicator, • semLm() to estimate linear models with an interval-censored dependent variable, • semLme() to estimate linear mixed models with an interval-censored dependent variable.
The functions are applied in order to estimate statistical indicators from interval-censored synthetic EU-SILC income data and to analyze interval-censored examination scores of students from London with linear and linear mixed regression models.
Further developments of the smicd package will include the possibility to estimate the bootstrapped standard errors in parallel computing environments. Additionally, it is planned to allow for the use of survey weights in the linear (mixed) regression models.