---
title: "Analysis of reporting behavior using generalized ordered probit models: an introduction to the *hopit* package"
author: "Maciej J. Dańko"
date: "29 III 2019"
output:
  pdf_document: 
    highlight: tango
    number_sections: no
    toc: no
fontsize: 11pt
geometry: margin=1.0in
abstract: The *hopit* package provides R functions to fit and analyze ordered response
  data in the context of reporting heterogeneity. In this vignette I describe the formulation
  and the fit of *hopit* models as well as functions that can be used to analyze heterogenity in reporting behavior. 
vignette: >
  %\VignetteIndexEntry{Introduction to hopit}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ../inst/REFERENCES_VIG.bib
biblio-style: apalike
link-citations: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
unlink('vignettes/vignette_cache', recursive = TRUE)
```

## 1. Introduction

*hopit* is an open source software library written in the R [@Rteam2019] and C++ [@RcppEigen2013; @Rcpp2011] programming languages. The *hopit* package provides versatile methods for fitting and analyzing ordered response data in the context of heterogeneity in self reporting behavior.

The ordered response data classify a measure of interest into ordered categories collected during a survey. For example, if the dependent variable is a happiness rating, then a respondent typically answers a question such as: “Taking all things together, would you say you are ... ?$"$; and then selects from response options such as: "very happy", "pretty happy", "not too happy", and "very unhappy" [@Liao2005]. Similarly, if interviewees are asked to evaluate their health in general (e.g., “Would you say your health is ... ?”), they may choose among several categories, such as "very good", "good", "fair", "bad", and "very bad" [@King2004; @Jurges2007; @Rebelo2014; @OKSUZYAN2019]. In political science, a respondent may be asked for an opinion about recent legislation (e.g. “Rate your feelings about the proposed legislation.") and asked to choose among categories like "strongly oppose", "mildly oppose", "indifferent", "mildly support", and "strongly support" [@GreeneHensher2010]. It is easy to imagine other multi-level ordinal variables that might be used during a survey and to which the methodology described below could be applied.

In practice, it is assumed that when responding to a survey question about their general happiness, health, feelings, attitudes or other status, participants are assessing their true value of this unobserved continuous variable, and project it onto the discrete scale provided. The thresholds that individuals use to categorize their true status by selecting a specific response option may be affected by the reference group chosen, their earlier life experiences, and cross-cultural differences in using scales. Thus, the responses of individuals may differ depending on their gender, age, cultural background, education, and personality traits; among other factors.

From the perspective of reporting behavior modeling, one of the main tasks researchers face is to compute this continuous estimate of the underlying, latent measures of individuals based on several specific characteristics of the responses considered (e.g., health variables or happiness variables), and to account for variations in reporting across socio-demographic and cultural groups. More specifically, to build a latent, underlying measure, a generalized hierarchical ordered threshold model is fitted that regresses the reported status/attitude/feeling on two sets of independent variables [@Boes2006; @Green2014]. When the dependent reported ordered variable is self-rated health status, then the first set of variables – i.e., health variables – assess specific aspects of individuals’ health, such as measures of chronic conditions, mobility, difficulties with a range of daily activities, grip strength, anthropometric characteristics, and lifestyle behaviors. Using the second set of independent variables (threshold variables), the model also adjusts for differences across socio-demographic and cultural groups, such as differences in cultural background, gender, age, and education [@King2004; @Jurges2007; @OKSUZYAN2019; but see @Rebelo2014].

Once the model is fitted, its estimates (latent measure and threshold coefficients) can be used to calculate the differences in reporting behavior among groups of people with different contextual characteristics through the calculation of differences between the expected and the reported ordinal response measures [@Jurges2007; @OKSUZYAN2019]. 

  **Table 1.** Glossary.
```{r, echo=FALSE, results='asis'}
library(pander)
dd <- c(
'Categorical response           ',	'$y$',	'Dependent variable obtained during the survey',	'Self-rated health, self-rated happiness',
'Latent measure                 ',	'$h$',	'Modeled continuous latent measure of the investigated response variable',	'Latent health, latent happiness',
'Latent index                   ',	'$H$',	'Standardized latent measure',	'Health index, happiness index',
'Latent variables               ',  '  ---',	 	'Variables used to model the latent measure',	'Health variables, happiness variables',
'Latent terms                   ',	'$X$',	'Terms of the design matrix used to model the latent measure', '',
'Latent coefficients            ',	'$\\beta$',	'Coefficients corresponding to each latent term','',
'Standardized coefficient           ',	'$D$',	'Standardized value of a coefficient', 'Disability weights',
'Thresholds                     ',  '$\\alpha$',	'Thresholds used to group the latent measure ', 	'Cut-points',
'Threshold variables            ',  '  ---',		'Variables used to model the thresholds	','Socio-demographic, cultural, contextual variables',
'Threshold terms                ',	'$Y$',	'Terms of the design matrix used to model the latent measure','',
'Threshold coefficients         ',	'$\\gamma$, $\\lambda$',	'Coefficients corresponding to each threshold term','')

dd <- matrix(dd, 11, 4, byrow=TRUE)
colnames(dd) <- c('Term',	'Symbol',	'Definition',	'Exemplary case specific synonyms')
dd <- as.data.frame(dd)
pander(dd,style='multiline',split.tables=Inf,justify='lccc',split.cells=c(Inf,3,25,25))
```

## 2. Generalized (hierarchical) ordered threshold model 

Ordered threshold models are used to fit ordered categorical dependent variables. The generalized ordered threshold models [@Terza1985; @Boes2006; @Green2014] are an extension of ordered threshold models [@McKelvey1975]. Wherease in the latter models the thresholds are constant, in the generalized models the thresholds are allowed to be dependent on covariates. @GreeneHensher2010 and @Green2014 pointed out that for a model to make a sense, the thresholds must also be ordered. This observation motivated Greene and coauthors to call these models *HOPIT*, which stands for hierarchical ordered probit models.

In the self-rated health example, the response variable is self-rated health and the latent measure $h_i$ can depend on different health conditions and diseases (health variables). These variables are modeled with the parallel regression assumption. According to this assumption, the coefficients that describe the relationship between the lowest response category and all of the higher response categories, are the same as the coefficients that describe the relationship between another (e.g. adjacent) the lowest response category and the remaining higher response categories. In the considered case $h_i$ is modeled as a linear function of the design matrix of health variables $X$ and its corresponding coefficients $\beta$:
\begin{equation}
\label{eq:1}
h_{i} = \sum_{k=1}^K \beta_kX_{i,k} = X'\beta
\end{equation}
where index $i \in 1...N$ is a number of cases (e.g. respondents), $X$ is in the form of a design matrix, and $K$ is number of columns in $X$.
As described above, the categorization (response mechanism) of the latent measure $h_i$ is modeled in terms of thresholds $\alpha_{i,j}$, while assuming that the lower order thresholds are never greater than the higher order thresholds (hierarchical assumption):

\begin{equation}
\label{eq:2}
\begin{cases}
y_i = 1 ~\Leftrightarrow~ \alpha_{i,0} \leq h_i < \alpha_{i,1}\\
y_i = 2 ~\Leftrightarrow~ \alpha_{i,1} \leq h_i < \alpha_{i,2}\\
\cdots\\
y_i = j~ \Leftrightarrow~ \alpha_{i,j-1} \leq h_i < \alpha_{i,j}\\
\cdots\\
y_i = J~ \Leftrightarrow~ \alpha_{i,J-1} \leq h_i < \alpha_{i,J}\\
\end{cases}
\end{equation}

The thresholds (cut points, $\alpha$) are modeled using threshold variables coded as design matrix $Y$, their coefficients $\gamma$, and intercepts $\lambda$. It is assumed that they model the contextual characteristics of the respondent (e.g. country, gender, and age). The threshold variables are modeled without applying the parallel regression assumption, thus each threshold is modeled by a variable independently [@Boes2006; @Green2014]. 

Different parametrizations of thresholds exist [@Green2014; @Rebelo2014; @Jurges2007]. In the package, parametrization of @King2004 and @Jurges2007 is used, which assumes that: 

\begin{equation}
\label{eq:3}
\alpha_{i,~j} =\begin{cases} -\infty& for~j=0 \\
  \lambda_{1} + \sum_{m=1}^{M} \gamma_{1,m} Y_{i,m} &for~j=1\\ 
 \alpha_{i,~j-1} +exp(\lambda_{j}+\sum_{m=1}^M \gamma_{j,m} Y_{i,m})&for~J-1 \ge j\ge2\\
 \infty& for~j=J
\end{cases}
\end{equation}

The condition $y_i = j~ \Leftrightarrow~ \alpha_{j-1,i} \leq h_i < \alpha_{j,i}$ can be easily expressed in terms of the probability, which leads to:
\begin{equation}
\label{eq:4}
P(y_i = j) = P(\alpha_{j-1,i} \leq h_i < \alpha_{j,i}),
\end{equation}
hence
\begin{equation}
\label{eq:5}
P(y_i = j) = \Phi(\alpha_{i,~j}-h_{i})-\Phi(\alpha_{i,~j-1}-h_{i}),
\end{equation}

where $\Phi$ is a distribution function (cdf, cumulative density function). For example, for probit regression it is standard normal cdf $\Phi(x)=\frac{1}{2}+\frac{1}{2}*erf \Big(\frac{x}{\sqrt 2}\Big)$ whereas for logit regression it takes the form $\Phi(x)=\frac{1}{1+e^{-x}}$. In reporting behavior analyses the probit model is typically chosen. This model simply assumes that $h_i$ is affected by a random noise $\epsilon_i$ having standard normal distribution $\epsilon_i\sim \mathcal{N}(0,1)$. 

Using all definitions presented above, the log likelihood function can be constructed as:
\begin{equation}
\label{eq:6}
\ln L = \sum_{i=1}^N \sum_{j=1}^J z_{i,~j} \ln\Big[\Phi(\alpha_{i,~j}-h_{i})-\Phi(\alpha_{i,~j-1}-h_{i})\Big],
\end{equation}
where $z_{i,j}$ is an indicator function defined as:
\begin{equation}
\label{eq:7}
z_{i,~j} =\begin{cases} 0&for~ y_i=j\\ 1&for~y_i\ne j\end{cases}
\end{equation}

## 3. Analysis of reporting heterogeneity

The model estimates are used to determine reporting behavior; i.e., how the continuous latent measure is projected onto the categorical response. In practice, this is done by comparing actual categorical ordered responses with theoretical responses that are adjusted for heterogeneity in reporting behaviors, and are, therefore, more comparable across individuals.

One of the first steps of the analysis is the standardization of the latent measure to obtain the latent index $H_i$.

\begin{equation}
\label{eq:8}
H_i = 1-\frac{h_i-\displaystyle\min_i h_i}{\displaystyle\max_i h_i-\displaystyle\min_i h_i}
\end{equation}

In the self-rated health example $H_i$ is a proxy for the true underlying health of an individual, and varies from 0, representing the (model-based) worst health state in the sample, to 1, representing the (model-based) best health state in the sample. 

The predicted latent measure $h_i$ obtained from the model is also used to standardize the latent variable coefficients. In the self-rated health example the standardized coefficients are called disability weights $D_k$ [@Jurges2007; @OKSUZYAN2019] and are calculated for each health variable to provide information about the impact of a specific health measure on the latent index $H_i$. The disability weight for a health variable is equal to the ratio of the corresponding health coefficient and the difference between the lowest and the highest values of predicted latent health. In other words, the disability weight reduces $H_i$ by some given amount or percentage (i.e. the $H_i$ of every individual is reduced by the same amount if the person had a heart attack or other heart problems)[@Jurges2007; @OKSUZYAN2019].

\begin{equation}
\label{eq:9}
D_k= \frac{\beta_k}{\displaystyle\max_i h_i-\displaystyle\min_i h_i}
\end{equation}

While the latent index $H_i$ is intended to reflect the underlying health, happiness or other status across individuals, the standardized coefficients $D_k$ (e.g. disability weights), are computed for an average individual in the study population. The relationship between $H_i$ and $D_k$ follows the equation:

\begin{equation}
\label{eq:10}
H_i = C-\sum_{k=1}^K D_kX_{i,k}, ~~~\text{where}~C=\frac{\displaystyle\max_i h_i}{\displaystyle\max_i h_i-\displaystyle\min_i h_i}
\end{equation}

Reporting behavior analysis is based on the reclassification of individuals into new response categories. There are two methods of reclassification: (1) @Jurges2007 percentile method [see also @Rebelo2014 and @OKSUZYAN2019] and (2) reclassification based on estimated thresholds.  

In the first method, the classification is based on the calculated latent index $H_i$ and is thus adjusted for inter-individual differences in reporting behavior. Jürges' percentile method is based on the original distribution of the categorical response variable. First for each category $j$ an empirical distribution function is constructed.

\begin{equation}
\label{eq:11}
\hat{F}(j) = \frac{1}{N}\sum_{i=1}^N \textbf{1}_{y_i} \leq j
\end{equation}

where $\textbf{1}$ is an indicator function taking 1 if the condition is true, or is 0 otherwise. The calculated cumulative frequencies of the latent index $H_i$ are used as percentiles (cut-points), so that each individual $i$ can be reclassified into new response categories.

In the second method the reclassification is based on eq. (2), so that each individual has her/his own, model-derived cut-points.

## 4. Installing and loading the package

The newest version of the package is always available from GitHub. It can be installed using the *devtools* package
```{r, echo=TRUE, eval=FALSE}
library(devtools)
install_github("maciejdanko/hopit")
```

```{r, echo=TRUE, eval=FALSE}
library(hopit)
```

```{r, echo=FALSE,  results='hide', eval=TRUE, include=FALSE}
g <- capture.output(library(hopit))
```
In the examples presented below I use *healthsurvey*, which is a completely artificial data set yhat is simulated using the distributions of some major health and socio-demographic characteristics. The distributions and the data structure is roughly based on the WAVE1 SHARE database (DOIs: 10.6103/SHARE.w1.600) see @Borsch2013 for technical details. See also the acknowledgements.

```{r, echo=TRUE, cache=TRUE}
# load *healthsurvey* dataset
data(healthsurvey)

# horizontal view of the dataset (omitting ID)
print(t(healthsurvey[1:6,-1]), quote=FALSE, na.print='NA', right=TRUE)
```

The first variable on the list (*health*) is the categorical self-reported health status. This variable is followed by 11 determinants of health, which includes information on the presence of chronic diseases and other health conditions. The *sex*, *ageclass*, *education*, and *country* variables describe the contextual characteristics of individuals. The last group of variables (*csw*, *psu*, and *ssu*) describes the survey design.

## 5. Fitting the model using the *hopit*() function

The generalized ordered probit model can be fitted using the *hopit*() function. The function takes two kinds of formulas: (1) *latent.formula*, which models the impact of the latent variables on categorical health; and (2) *thresh.formula, which models the thresholds.   

```{r, echo=TRUE, cache=TRUE} 
# first determine the order of the dependent variable
levels(healthsurvey$health)

# the order is decreasing (from the best to the worst health state)
# so we set: decreasing.levels = TRUE
model1<- hopit(latent.formula = health ~ hypertension + high_cholesterol + 
                             heart_attack_or_stroke + poor_mobility + very_poor_grip + 
                             depression + respiratory_problems + 
                             IADL_problems + obese + diabetes + other_diseases, 
               thresh.formula = ~ sex + ageclass,
               decreasing.levels = TRUE,
               control=list(trace=FALSE),
               data = healthsurvey)

summary(model1)
```

*model1* contains 11 dichotomous health variables and two threshold variables. The fitted coefficient includes $\beta$s ("latent.params", first 11 coefficients in the summary), $\lambda$s ("thresh.lambda", threshold intercepts, "(L)" prefix in the summary), and $\gamma$s ("thresh.gamma", parameters related to threshold covariates, "(G)" prefix in the summary). The model coefficients can be accessed by the *coef*() function.

```{r, echo=TRUE, cache=TRUE}
# extract parameters in the form of a list
cm1 <- coef(model1, aslist = TRUE)

# names of the returned coefficients
names(cm1)

# extracting the latent health coefficients
cm1$latent.params
```

*model1* can be further extended by adding the country of origin to the threshold formula to control for cultural differences.

```{r, echo=TRUE, cache=TRUE}
model2<- hopit(latent.formula = health ~ hypertension + high_cholesterol + 
                      heart_attack_or_stroke + poor_mobility + 
                      very_poor_grip + depression + respiratory_problems + 
                      IADL_problems + obese + diabetes + other_diseases, 
               thresh.formula = ~ sex + ageclass + country,
               decreasing.levels = TRUE,
               control=list(trace=FALSE),
               data = healthsurvey)
```

The fit of the two models can be compared using the AIC() function:

```{r, echo=TRUE, cache=TRUE}
AIC(model2, model1)
```

or using the Likelihood Ratio Test (LRT) because the models are nested:

```{r, echo=TRUE, cache=TRUE}
anova(model2, model1)
```

Both *latent.formula* and *thresh.formula* allow the user to specify the interactions, like the interaction between gender (*sex*) and age (*ageclass*): 

```{r, echo=TRUE, cache=TRUE}
model3<- hopit(latent.formula = health ~ hypertension + high_cholesterol + 
                      heart_attack_or_stroke + poor_mobility + 
                      very_poor_grip + depression + respiratory_problems + 
                      IADL_problems + obese + diabetes + other_diseases, 
               thresh.formula = ~ sex * ageclass + country,
               decreasing.levels = TRUE,
               control=list(trace=FALSE),
               data = healthsurvey)

print(anova(model3,model2), short=TRUE)

```

The interactions between latent and threshold variables can also be modeled. Depending on how an interaction is interpreted, it can be added to either the latent or the threshold formula:

```{r, echo=TRUE, cache=TRUE}
model4<- hopit(latent.formula = health ~ hypertension + high_cholesterol + 
                      heart_attack_or_stroke + poor_mobility + 
                      very_poor_grip + depression + respiratory_problems + 
                      IADL_problems + obese + diabetes + other_diseases +
                      sex : respiratory_problems, 
               thresh.formula = ~ sex * ageclass + country + sex : depression,
               decreasing.levels = TRUE,
               control=list(trace=FALSE),
               data = healthsurvey)

print(anova(model3,model4), short=TRUE)

```

The *hopit*() function also has an option to include the survey design using the *survey* package. 
The example below fits a model using a simple two-level cluster sampling design.

```{r, echo=TRUE, cache=TRUE}
design <- svydesign(ids = ~ country + psu, weights = healthsurvey$csw, 
                    data = healthsurvey)

model2s<- hopit(latent.formula = health ~ hypertension + high_cholesterol + 
                       heart_attack_or_stroke + poor_mobility + 
                       very_poor_grip + depression + respiratory_problems + 
                       IADL_problems + obese + diabetes + other_diseases, 
                thresh.formula = ~ sex + ageclass + country,
                decreasing.levels = TRUE,
                design = design,
                control=list(trace=FALSE),
                data = healthsurvey)
```

Generally, ignoring the survey design can lead to biased results. In the example presented here, the design has little importance, which is seen by comparing the coefficients of the latent variable for the two models:

```{r, echo=TRUE, cache=TRUE}
cbind('No survey design'=coef(model2,aslist=TRUE)$latent.par,
      'Has survey design'=coef(model2s,aslist=TRUE)$latent.par)
```

The accuracy of the model fit can be assessed using the *profile*() function, which calculates and plots the profile of the log likelihood function around fitted coefficient values.

```{r, echo=TRUE, cache=TRUE}
profile(model3)
```

## 6. Analyses of reporting heterogeneity using the *hopit* package

Let us now look at the latent health variables in *model3*. 

```{r, echo=TRUE, cache=TRUE}
model3$coef.ls$latent.params
```

We can standardize the coefficients of these variables using Jürges' approach [@Jurges2007] to obtain disability weights. The standardization can be done using the *standardizeCoef*() function.

```{r, echo=TRUE, fig.height = 5, fig.width = 5, fig.align = "center", cache=TRUE}
# A function that modifies the coefficient names.  
txtfun <- function(x) gsub('_',' ',substr(x,1,nchar(x)-3))

# Calculate and plot the disability weights
sc <- standardizeCoef(model3, plotf = TRUE, namesf = txtfun)
sc
```

The *namesf* argument is a function or a character vector that is used to rename the coefficients. Here, it removes the last 3 letters ("yes"), which is a reference level for each variable, and replaces "_" with spaces in variable names.

The latent index is simply calculated using the *latentindex*() function.

```{r, echo=TRUE, fig.height = 4, fig.width = 5, fig.align = "center", cache=TRUE}
hi <- latentIndex(model3, plotf = TRUE, response = "data", 
                  ylab = 'Health index', col='deepskyblue3')
```

The boxplot above shows the reported health status versus the health index. It is also possible to plot the expected categorical health status on the Y axis calculated according to eq. (2): 

```{r, echo=TRUE, fig.height = 4, fig.width = 5, fig.align = "center", cache=TRUE, eval=FALSE}
hi <- latentIndex(model3, plotf = TRUE, response = "fitted", 
                  ylab = 'Health index', col='deepskyblue3')
```

or according to @Jurges2007 method:

```{r, echo=TRUE, fig.height = 4, fig.width = 5, fig.align = "center", cache=TRUE, eval=FALSE}
hi <- latentIndex(model3, plotf = TRUE, response = "Jurges", 
                  ylab = 'Health index', col='deepskyblue3')
```

The main aim of reporting heterogeneity analyses is to determine the cut-points used to calculate the adjusted health status for each individual. The calculation and the plotting of cut-points is realized using the *getCutPoints*() function.

```{r, echo=TRUE, fig.height = 3.8, fig.width = 4.6, fig.align = "center", cache=TRUE}
z <- getCutPoints(model=model3)

# Health index cut-points
z$cutpoints

# Adjusted health levels for individuals (using the Jürges method)
rev(table(z$adjusted.levels))

# Original health levels for individuals
table(model3$y_i)

# Adjusted health levels for individuals (using estimated model thresholds)
table(model3$Ey_i)
```

The analysis of health levels is performed using the *getLevels*() function

```{r, echo=TRUE, cache=TRUE, fig.height = 4, fig.width = 6, fig.align = "center", cache=TRUE}
# Health levels for age and gender, and pooled country of origin.
hl <- getLevels(model = model3, formula = ~ sex + ageclass, data = healthsurvey, 
                      sep=' ', plotf=TRUE)
```

The differences between the original and the adjusted frequencies can be calculated directly using the *getLevels* output:

```{r, echo=TRUE}
round(100*(hl$original - hl$adjusted),2) # in (%)

```

## 7. Bootstrapping Confidence Intervals

The package offers functions for calculating the confidence intervals for any measure derived from the model. In the example below, we calculate the confidence intervals of the difference between the original and the adjusted frequencies of bad health. The adjusted frequencies are determined by the presence of "*Poor*" or "*Fair*" self-rated health categories.

```{r, echo=TRUE, fig.height = 5.6, fig.width = 4.7, fig.align = "center", cache=TRUE}
# the function to be bootstraped
diff_BadHealth <- function(model, data) {
  hl <- getLevels(model = model, formula = ~ sex + ageclass, data = data, 
                  sep = ' ', plotf = FALSE)
  hl$original[,1] + hl$original[,2] - hl$adjusted[,1]- hl$adjusted[,2]
}

# estimate the difference
est.org <- diff_BadHealth(model = model3, data = healthsurvey)

# perform the bootstrap
B <- boot_hopit(model = model3, data = healthsurvey, 
                func = diff_BadHealth, nboot = 100)

# calculate lower and upper bounds using the percentile method
est.CI <- percentile_CI(B)

# plot the difference and its (asymmetrical) confidence intervals
pmar <- par('mar'); par(mar = c(9.5,pmar[2:4]))
m <- max(abs(est.CI))
pos <- barplot(est.org, names.arg = names(est.org), las = 3, ylab = 'Original - Adjusted', 
               ylim=c(-m, m), density = 20, angle = c(45, -45), col = c('blue', 'orange'))
for (k in seq_along(pos)) lines(c(pos[k,1],pos[k,1]), est.CI[,k], lwd = 2, col = 2)
abline(h = 0); box(); par(mar = pmar)

```

The results show that men tend to over-report bad health at ages (50,60] and (50,70], whereas women tend to over-report bad health at ages [70,80); and that both sexes at ages (80, 120] tend to under-report bad health.

## 8. Acknowledgements

I thank Anna Oksuzyan, Christian Dudel, Marius Pascariu, Laszlo Nemeth, and Oskar Burger for their help and suggestions on the package and the vignette. I also thank the Max-Planck Institute for Demographic Research for all their support.

The artificially generated data used in the package examples is based on the distributions of some major socio-demographic and health related characteristics of the WAVE 1 SHARE database. None of the records represent any part of the true data. 

The SHARE data collection has been primarily funded by the European Commission through FP5 (QLK6-CT-2001-00360), FP6 (SHARE-I3: RII-CT-2006-062193, COMPARE: CIT5-CT-2005-028857, SHARELIFE: CIT4-CT-2006-028812) and FP7 (SHARE-PREP: N°211909, SHARE-LEAP: N°227822, SHARE M4: N°261982). Additional funding from the German Ministry of Education and Research, the Max Planck Society for the Advancement of Science, the U.S. National Institute on Aging (U01_AG09740-13S2, P01_AG005842, P01_AG08291, P30_AG12815, R21_AG025169, Y1-AG-4553-01, IAG_BSR06-11, OGHA_04-064, HHSN271201300071C) and from various national funding sources is gratefully acknowledged (see www.share-project.org).

## 9. References

