longtable

# A Journey from Wild to Textbook Data to Reproducibly Refresh the Wages Data
from the National Longitudinal Survey of Youth Database

Dewi Amaliah
Department of Econometrics and Business Statistics
Monash University
Clayton, VIC 3800
<EMAIL_ADDRESS>
Dianne Cook
Department of Econometrics and Business Statistics
Monash University
Clayton, VIC 3800
<EMAIL_ADDRESS>
Emi Tanaka
Department of Econometrics and Business Statistics
Monash University
Clayton, VIC 3800
<EMAIL_ADDRESS>
Kate Hyde
Department of Econometrics and Business Statistics
Monash University
Clayton, VIC 3800
<EMAIL_ADDRESS>
Nicholas Tierney
Telethon Kids Institute
Nedlands, WA 6009
<EMAIL_ADDRESS>

###### Abstract

Textbook data is essential for teaching statistics and data science methods
because they are clean, allowing the instructor to focus on methodology.
Ideally textbook data sets are refreshed regularly, especially when they are
subsets taken from an on-going data collection. It is also important to use
contemporary data for teaching, to imbue the sense that the methodology is
relevant today. This paper describes the trials and tribulations of refreshing
a textbook data set on wages, extracted from the National Longitudinal Survey
of Youth (NLSY79) in the early 1990s. The data is useful for teaching modeling
and exploratory analysis of longitudinal data. Subsets of NLSY79, including
the wages data, can be found in supplementary files from numerous textbooks
and research articles. The NLSY79 database has been continuously updated
through to 2018, so new records are available. Here we describe our journey to
refresh the wages data, and document the process so that the data can be
regularly updated into the future. Our journey was difficult because the steps
and decisions taken to get from the raw data to the wages textbook subset have
not been clearly articulated. We have been diligent to provide a reproducible
workflow for others to follow, which also hopefully inspires more attempts at
refreshing data for teaching. Three new data sets and the code to produce them
are provided in the open source R package called yowie.

_K_ eywords Data cleaning; Data tidying; Reproducible workflow; Longitudinal
data; NLSY79; Initial data analysis;

##  1 Introduction

Statistics and data science education relies on cleaned and simplified data,
suitably called textbook data, for clear examples about how to apply different
techniques. An example of this is the wages data made public by Singer and
Willett (2003) in their book, “Applied longitudinal data analysis,” which can
be used to teach generalized linear models, and additionally hierarchical,
mixed effects, and multilevel models. The data records hourly wages of a
sample of high school dropouts from 1979-1994, along with the demographic
variables, such as education and race, taken from the National Longitudinal
Survey of Youth (NLSY79) (Bureau of Labor Statistics, U.S. Department of Labor
2021a).

The story from modeling the data (and as reported by Singer and Willett) is
that wages increase with the length of time in the workforce, a higher level
of education leads to higher wages, and that race makes a difference, on
average. However, an exploratory analysis reveals that the individual
experience varies a great deal from the overall average. Some individuals
experience a decline in wages the longer they are in the workforce, and many
experience volatility in their wages. It is for these reasons that the wages
data was used to illustrate exploratory longitudinal data analysis in Ilk
(2004), and was further developed into a case study for use in the teaching of
exploratory data analysis at Iowa State University.

This disparity between the overall trend and the individual heterogeneity is
what makes this data interesting. Textbook data sets have longevity if they
have unresolved elements. The iris data (Anderson 1935) is a prime example. It
has withstood the test of time because the three species cannot be perfectly
classified, and so it continues to challenge researchers and instructors to do
better in the analysis111It should be noted that the iris data is best
replaced today with the penguins data (Horst, Hill, and Gorman 2020), which
has similar qualities, is new, and does not suffer from a connection with
eugenics (Stodel 2020). We argue that the wages data is in this class of
textbook data, too, because it presents a challenge for longitudinal data
analysis: how can we better summarize and explain the individual experience?

For statistics, and data science by association, it is increasingly important
to reach the individual. One might think of this as a divergence of purpose –
statistics for public policy – or statistics for the public. The two are not
the same. As the world becomes more electronically connected, combating
misinformation and mitigating conspiracy theories require statistics that
address the individual. For example, with the wages data, the overall trend,
across demographics, is a steady increase in wages, but the patterns among
individuals is heterogeneous. Some individuals see a decline in their wages,
some an increase, some have dramatic ups and downs from year to year, and the
heterogeneity of the patterns are evident across demographics. Subject-
specific variation explains most of the variation in wages, and even though
the overall trend is statistically significant, it is weak. The message for
public policy is that demographic profile is related to different wage
patterns on average and some structural social change is desirable. However,
discussing the overall trend with an individual is misleading – more than
likely their pattern is completely different. It could even be depressing for
them to compare themselves to the overall average, especially if the
individual’s wages have a declining trend or high volatility. It would be more
useful for an individual to know what percentage of the population have
profiles like theirs. We argue that different summaries are needed to
communicate with the public, and ones where they might be able to individually
act on to change their own situation. We occasionally hear this voiced in the
public media, also, for example, an article published in the Sydney Morning
Herald argues there is no average Australian (Moncrief 2015). Thus the wages
data provides educators with a challenge, from a methodological and a
philosophical perspective, and wage experience is a topic of interest to many.

As a textbook data set, though, the wages data is outdated. The most recent
year in the data is 1994, 9 years prior to when Singer and Willett (2003) was
published. Teachers of statistics need contemporary data sets to show how
techniques are relevant for today’s students. Using tired old textbook data
sets can imbue a misconception that the field is not current. The wages data
is extracted from NSLY79, one of the best examples of open data (see details
at Open Knowledge Foundation 2021), which is constantly being updated. It
should be possible to continuously refresh the textbook data from the data
repository. This paper describes our (non-glamorous) journey from open wild
data to textbook data.

This paper demonstrates the steps of cleaning data, including subjective
decisions made on dealing with anomalies, and documents the process, as
recommended by Huebner, Vach, and Cessie (2016). They emphasize that making
the data cleaning process accountable and transparent is imperative and
essential for the integrity of downstream statistical analyses and model
building. Clean data often then goes through an “initial data analysis” (IDA)
(Chatfield 1985), where one would summarize and scrutinize the data,
especially to check if the data is consistent with assumptions required for
modeling. This stage is related to exploratory data analysis (EDA), coined by
Tukey (1977) with a focus on learning from data. EDA can be considered to
encompass IDA. In practice, the three stages of cleaning, summarizing, and
exploring are cyclical, one often needs to do more cleaning after
scrutinizing. Dasu and Johnson (2003) say that data cleaning and exploration
is a difficult task and typically consumes a large percentage of the time
spent in analyzing data.

Our approach to cleaning builds heavily on the tidyverse approach (Wickham et
al. 2019). The data is first organized into “tidy data” (Wickham 2014) and
then further wrangled using step-wise piping with a split-apply-combine
strategy for mutating new variables (Wickham 2011). Tidy data shouldn’t be
confused with “tame data,” which Kim, Ismay, and Chunn (2018) coined to refer
to textbook data sets suitable for teaching, particularly teaching statistics.
The resulting (tame) data is provided in a new R package called yowie, which
includes the code so that the process is reproducible and could be used to
further refresh the data as new records are made available in the NLSY79
database.

This paper is structured in the following way. Section 2 describes the NLSY79
data source. Section 3 presents the steps of cleaning the data, including
getting and tidying the data from the NLSY79 and IDA to find and repair
anomalies. Our final subset is compared to the old textbook subset in Section
4. Finally, Section 5 summarizes the contribution and makes recommendations
for the NLSY79 data curators.

##  2 The NLSY79

Singer and Willett (2003) used the wages and other variables of high school
dropouts from the NLSY79 data as an example data set to illustrate
longitudinal data modeling of wages on workforce experience, with covariates
education and race. This data has been playing an important role in research
in various disciplines, including but not limited to economics, sociology,
education, public policy, and public health for more than a quarter of the
century (Pergamit et al. 2001). In addition, this is considered a carefully
designed longitudinal survey with high retention rates, making it suitable for
life course research (Pergamit et al. 2001; Cooksey 2017). According to
Cooksey (2017), thousands of articles and hundreds of book chapters and
monographs have utilized this data. Moreover, the NLSY79 is considered the
most widely used and most important cohort in the survey data (Pergamit et al.
2001).

Our aim is to refresh the wages textbook data and append it with data from
1994 through to the latest data reported in 2018, a purpose consistent with
Grimshaw (2015)’s statistics education goal of embracing authentic data
experiences. Here, we investigate the process of getting from the raw NLSY79
data to a textbook data set as similar as possible to that provided by Singer
and Willett (2003). We should also note that race is a variable in the
original data set, and for compatibility, it is also provided with the
refreshed data for the _purposes of studying racism, not race_ (Fullilove
1998). There are a number of data sets provided by Singer and Willett (2003),
but we focus only on wages data because it has captivated our attention for a
number of years in our own teaching of longitudinal data analysis.

###  2.1 Database

The NLSY79 is a longitudinal survey administered by the U.S Bureau of Labor
Statistics that follows the lives of a sample of American youth born between
1957-1964 (Bureau of Labor Statistics, U.S. Department of Labor 2021a). The
cohort originally included 12,686 respondents aged 14-22 when first
interviewed in 1979. For a variety of reasons, some structural, the number of
respondents dropped to 9,964 after 1990. The surveys were conducted annually
from 1979 to 1994 and biennially thereafter. Data are currently available from
Round 1 (1979 survey year) to Round 28 (2018 survey year).

Although the main focus area of the NLSY is labor and employment, the NLSY
also covers several other topics, including education, training, achievement,
household, geography, dating, marriage, cohabitation, sexual activity,
pregnancy, fertility, children, income, assets, health, attitudes and
expectations, crime, and substance use.

There are two ways to conduct the interview of the NLSY79, which are face-to-
face or telephone interviews. In recent survey years, more than 90 percent of
respondents were interviewed by telephone (Cooksey 2017).

###  2.2 Target data

The NLSY79 data used in Singer and Willett (2003) contains the longitudinal
records of male high school dropouts who first participated in the study at
age 14-17 years from 1979 through to 1994. This dataset contains several
variables as follows:

  1. 1.

ID: the respondents’ ID.

  2. 2.

EXPER: stands for experience, the temporal scale, i.e., the length of time
(years) in the workforce, starting on the respondents’ first day at work.

  3. 3.

LNW: natural logarithm of wages, adjusted with 1990’s inflation rate.

  4. 4.

BLACK: binary variable, 1 indicates Black and 0 otherwise.

  5. 5.

HISPANIC: binary variable, 1 indicates Hispanic and 0 otherwise.

  6. 6.

HGC: the highest grade completed.

  7. 7.

UERATE: the unemployment rate of the year of the survey. When missing, the
variable is set to be 7.875 (the average rate).

We refresh this data by re-creating the full data with records from survey
years 1979 through to 2018 (the most recent year published). We also modify
some variables. For example, we use a single categorical race variable instead
of the two binary race variables, for reasons detailed below. We also include
additional variables, some for the purpose of providing more options for data
exploration in teaching examples: year of the survey, age of individual in
1979, whether the individual completed high school with a diploma or with a
graduate equivalency degree (GED), the highest grade completed in the
corresponding year of survey, the number of jobs the individual had in the
corresponding year of survey, the total number of hours the individual usually
works per week, the year when the individual started to work, and the number
of years the individual worked. We do not attempt to refresh the unemployment
rate variable.

We aim to create three datasets as follows:

  1. 1.

The wages data of the whole NLSY79 cohort, including females,

  2. 2.

A separate table of the demographic data of the whole NLSY79 cohort, and

  3. 3.

The subset of wages data in (1), which is the wages of high school dropouts’
as a refreshed version of Singer and Willett (2003)’s data.

##  3 Data cleaning

In the context of official statistics, M. van der Loo and de Jonge (2018)
describe the “statistical value chain,” which includes various production
stages of the data cleaning process as raw data (initial data as delivered
originally), input data (data organized with correct type and identified
variables), and valid data (data that has been cleaned and more accurately
represents the intent of variables). What we have colorfully named wild data
can be considered to be raw data, and valid data could be considered to be
textbook data in the above statistical value chain. In this section, we
outline the steps to download the raw data (Section 3.1) and then tidy the raw
data into input data, specifically for the demographic variables (Section
3.1.1) and the employment variables (Section 3.1.2), so that the resulting
input data can be used downstream for validating the data as described in
Section 3.3 and Section 3.4.

###  3.1 Getting the data

The NLSY79 data contains a large number of variables, but for our purposes,
the scope required is limited to demographic profiles, wages data, and work
experience. More specifically, we went to the NLSY79 database website at
https://www.nlsinfo.org/content/cohorts/nlsy79/get-data, clicked on the direct
link to NSLY79 data, and navigated as described in Figure 1.

Navigating the data source NLSY79
(https://www.nlsinfo.org/investigator/pages/search?s=NLSY79) The CASEID will
be always be selected. The 3 recommended demographic variable (sample ID, race
and sex) were selected. For the remaining variables, we went to the “Variable
Search" tab and select variables as follows $\triangleright$ Education,
Training and Achievement Scores $\triangleright$ Education $\triangleright$
Summary measures $\triangleright$ All schools $\triangleright$ By year
$\triangleright$ Highest grade completed All 80 variables in Highest grade
completed were selected. $\triangleright$ Dates of diploma or degree All
variables named Q3-8A were selected. $\triangleright$ Employment
$\triangleright$ Summary measures $\triangleright$ By job $\triangleright$
Hours worked All 447 primary variables in Hours worked were selected.
$\triangleright$ Hourly wages All 156 variables in Hourly wages were selected.
$\triangleright$ Summary measures $\triangleright$ Since date of last
interview $\triangleright$ Weeks worked All 28 variables in Weeks worked were
selected. $\triangleright$ Employer Roster $\triangleright$ Job dates
$\triangleright$ Original start date Only selected the start date (Year) for
the first job (E00101.02) $\triangleright$ Household, Geography & Demographics
$\triangleright$ Demographics $\triangleright$ Basic demographics
$\triangleright$ Date of birth All 4 variables in Date of birth were selected.
To download all 742 variables selected, we then navigate to the tab “Save /
Download" then select the tab “Advanced Download". We select the R Source code
and comma-delimited datafile of selected variables with Reference Number as
column headers. We name the filename “NLSY79" and press the download button.
There are also options to get control or dictionary files for SAS, SPSS and
STATA. Figure 1: Documented steps taken to select variables of interest and
download the raw data.

The downloaded data set comes as a zip file, containing the following set of
files:

  * •

NLSY79.csv: comma separated value format of the response data,

  * •

NLSY79.dat: alternative text format of the response data,

  * •

NLSY79.NLSY79: tagset of variables that can be uploaded to the website to
recreate the data set, and

  * •

NLSY79.R: R script for reading the data into R and converting the variable
names and label into something more sensible.

We alter only the file path in NLSY79.R and run the script without any other
alteration. This results in the initial processing of the raw data into two
data sets, categories_qnames (where the observations are stored in
categorical/interval values) and new_data_qnames (the observations are stored
in integer form).

To get the data into tidy form (Wickham 2014), it needs to comply with three
rules: (i) each variable forms a column, (ii) each observation forms a row,
and (iii) each type of observational unit forms a table. The raw data,
new_data_qnames, violates these rules (i) and (ii) because information about
an individual’s multiple jobs over different years are in multiple columns.
The raw data consequently has a large number of columns (742 to be specific).
The values in the cell under the variables begin with HRP correspond to the
hourly wage in dollars. A glimpse of this data show the format problems:

    #> Rows: 12,686
    #> Columns: 11
    #> $ CASEID_1979 <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1~
    #> $ HRP1_1979   <int> 328, 385, 365, NA, 310, NA, NA, NA, 214, NA, 337~
    #> $ HRP2_1979   <int> NA, NA, NA, NA, 375, NA, NA, NA, NA, NA, 300, NA~
    #> $ HRP3_1979   <int> NA, NA, 275, NA, NA, NA, NA, NA, NA, NA, NA, NA,~
    #> $ HRP4_1979   <int> NA, NA, NA, NA, NA, 250, NA, NA, NA, NA, NA, NA,~
    #> $ HRP5_1979   <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ~
    #> $ HRP1_1980   <int> NA, 457, 397, NA, 333, 275, 300, 394, 200, 318, ~
    #> $ HRP2_1980   <int> NA, NA, 367, NA, NA, NA, NA, NA, NA, NA, NA, NA,~
    #> $ HRP3_1980   <int> NA, NA, 380, NA, NA, NA, 290, NA, NA, NA, NA, NA~
    #> $ HRP4_1980   <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ~
    #> $ HRP5_1980   <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ~

Thus, we re-arrange and wrangle the data into tidy data form, with columns
corresponding to individual ID, year, job number, wage in dollars, and the
demographic variables. This is done using the tidyverse suite of packages
(Wickham et al. 2019): tidyr (Wickham 2020) to pivot the data into long-form,
with dplyr (Wickham et al. 2020) and stringr (Wickham 2019) to create new
variables from the downloaded data from the database, and code levels of
factors by text wrangling. The long form of the data makes it possible to do
these data transformations efficiently, and it is an intermediate step towards
the final target data. The code for tidying the data are demonstrated at
https://numbats.github.io/yowie/articles/raw-to-input-data.html but also
described in the subsequent subsections.

####  3.1.1 Tidying demographic variables

Our final target data will include the demographic variables (variable names):
sex (sex), race (race), age (age_1979), highest grade completed reported in
each round of the survey (grade), highest grade completed ever reported (hgc),
highest grade completed in terms of years (as an integers), e.g., 9th grade =
9, 3rd-year college = 15, (hgc_i), highest grade completed in 1979 (hgc_1979)
and whether the graduate equivalency diploma is obtained (ged).

For sex and race, we have simplified the original names provided in the raw
data, SAMPLE_SEX_1979 and SAMPLE_RACE_78SCRN, respectively. However, it is
important to note that the use of “sex” is probably not correct. This
information is provided by the individual, and although only has two
categories, is more consistent with “gender” as defined in Heidari et al.
(2016). When using the dataset in the classroom, the educators might include
discourse on the use of terms “sex” and “gender”. Particularly, measuring
gender as a binary variable has the potential to fail to capture people who do
not identify themselves as either male or female or create measurement error
for people whose gender does not align with their sex classification. From a
statistical perspective, this can make it difficult to adjust survey
statistics to the population when gender is measured in different
categorization in different data sources. Kennedy et al. (2020) provides a
helpful discussion on these topics.

Similarly, there are issues for the race variable. The current US Federal
guidelines (Office of Management and Budget 1997) state that there should be
at least five categories for race, and an individual should be able to
identify as more than one race category. In addition, “Hispanic” or “Latino”
should be considered to be an ethnic group rather than a race. This level of
detail is not available in the database. There is one variable with only three
categories, and an individual can be a member of only one. Thus, the race
variable is inadequate by today’s standards, and educators should point to the
current Federal guidelines when using this variable, and keep in mind that the
purpose of any analysis is to study racism rather than race.

Furthermore, the object new_data_qnames, contains the variables Q1-3_A~Y_1979
and Q1-3_A~Y_1981, which records two versions of the birth year of the
respondent; this is also the case for the record of birth month (Q1-3_A~M_1979
and Q1-3_A~M_1981). The record contains two versions of birth year and birth
month as the survey recorded this in 1979 and 1981. We checked for consistency
between the two versions and found no discrepancy where the responses were
recorded in both 1979 and 1981. The age was then calculated using the birth
year.

The next step is processing the highest grade completed. There are several
ways to define this, and this should be reflected in the refreshed data to
give some flexibility for downstream analysis. The first one is the highest
grade ever completed is reported in the database and provided in the refreshed
data as hgc and hgc_i, encoded as a factor and integer variable, respectively.
For each individual, there is only one value of hgc and hgc_i. This variable
is obtained from new_data_qnames with the name HGC_EVER_XRND and stored in
year units (e.g., 10, 11, 12, 13, and so on), so the transformation is simply
giving a new column name as hgc_i and recoded as a factor to give hgc (e.g.,
10th grade, 11th grade, and so on). We decide to include both the factor and
integer variables to be used as examples in various classroom demonstrations
to explain different data types say. Besides this, it can also serve as a
demonstration of different type of data visualization ascribed by the data
type and for teaching the modeling of cross sectional data.

The second definition is the highest grade completed that is reported in each
round of the survey (provided in the refreshed data as grade). This value can
change over time, but it should only increase. This has been included in the
refreshed dataset to enable richer downstream analyses, for example, if one
would like to explore how the temporal changes in education level affect the
wages of individuals. Hence, different from hgc, this variable could be used
to demonstrate longitudinal and time series analysis in the classroom. This is
recorded in new_data_qnames as columns beginning with Q3-4 and HGC with year
as a suffix. In addition, it is also in the columns beginning with HGCREV
reflecting revised data. We chose to use these revised values because there
were fewer missing values indicating it had been more thoroughly checked. In
the survey years when the revised value were not reported, i.e., in 2012,
2014, 2016, and 2018, we used the unrevised version accordingly.

The third definition is the highest grade completed in 1979 (provided in the
refreshed data as hgc_1979), corresponding to the value from the first round
of the survey. This is calculated from the grade value in 1979. This best
reflects the grade when the individual left school and is included in order to
compare with the original data.

The next step is tidying to obtain ged. Along with hgc, ged is used to subset
the data to the high school dropouts as in the original data. The graduate
equivalency status is saved as a variable started with “Q3-8A” followed by the
year of the survey. Thus, we only separate the year and the GED status.
Although the GED status is asked in each round of the survey, we only retain
the latest status of one’s GED.

Finally, we get all of the demographic profiles of the NLSY79 cohort. We then
save this data as demog_nlsy79.

####  3.1.2 Tidying employment variables

Our target variables for the employment are to obtain respondent’s mean hourly
wage (wage), the number of jobs (njobs), the total hours of work per week
(hours) for each survey year, the year when individual starting to work
(stwork), the length of time (years) in workforce (yr_wforce), and work
experience measured as the number of years worked (exp). As the data only
reports up to 5 jobs for each respondent, the maximum number of jobs is capped
at 5.

From 1979 to 1987, new_data_qnames only contains one version of hours worked
per week for each job (in the variables with names starting with QES-52A).
From 1988 onward, we selected the total hours worked per week, including hours
working from home (QES-52D). However, in 1993, this variable was missing for
the first and last job, so we selected to use QES-52A instead. In addition,
2008 only had jobs 1-4 for the QES-52D variable, so we use only these.

The hourly wages are in the variables beginning with HRP in new_data_qnames.
As a respondent may have multiple jobs, the mean_hourly_wage is computed as a
weighted average of the hourly wage for each job with the number of hours
worked for each job as weights (provided that the information on the number of
hours is available); if the number of hours worked for any job is missing,
then the mean_hourly_wage is computed as a simple average of all available
hourly wages. Prior to computing the mean hourly wage, we undertook a number
of steps to treat extreme observations as described below:

  * •

If the hourly rate is recorded as 0, we set wage as missing, and

  * •

If the total hours of worked per week for the corresponding job is greater
than 84 hours, we set the wage and hours worked as missing.

The number of jobs (number_of_jobs) for each respondent per year is computed
from the number of non-missing values of hourly wage. In other words, even if
the information of hours worked exists for a particular observation, we do not
tally when the hourly wage is missing.

    #> # A tibble: 10 x 6
    #>       id  year mean_hourly_wage total_hours number_of_jobs is_wm
    #>    <int> <dbl>            <dbl>       <int>          <dbl> <lgl>
    #>  1     1  1979             3.28          38              1 FALSE
    #>  2     1  1981             3.61          NA              1 FALSE
    #>  3     2  1979             3.85          35              1 FALSE
    #>  4     2  1980             4.57          NA              1 FALSE
    #>  5     2  1981             5.14          NA              1 FALSE
    #>  6     2  1982             5.71          35              1 FALSE
    #>  7     2  1983             5.71          NA              1 FALSE
    #>  8     2  1984             5.14          NA              1 FALSE
    #>  9     2  1985             7.71          NA              1 FALSE
    #> 10     2  1986             7.69          NA              1 FALSE

For stwork variable, we only rename the column names from new_data_qnames.
This variable is then used to calculate the next variable, yr_wforce for each
survey round, which is the year of survey (year) minus the year of individual
started working (stwork). Finally, exp variable is derived from the number of
weeks worked since the last interview indicated as variable started with WKSWK
in new_data_qnames. To obtain the work experience since 1979, we calculate the
cumulative value. As the measurement unit is weeks, we convert this to
fractional year.

The employment and demographic variables are then joined. These data are
further filtered to the cohort who participated in at least three rounds in
the survey, the minimum observation recommended for longitudinal data in
Singer and Willett (2003). We also opt to restrict the minimum number of
observations so that the data can be used to demonstrate within-person
variation when it is used to teach longitudinal data. However, it is worth
noting that the original data does not restrict the number of observations for
each individual, i.e., there are individuals with only one and two
observations.

Finally, we save the resultant wage data on this cohort as wages. Note that we
save grade variable in this dataset instead of demog_nlsy79 dataset because it
is a longitudinal variable, while demog_nlsy79 is a cross-sectional dataset
reflecting the state of individuals corresponding to the most recent round of
the survey.

###  3.2 Calculated variables: work experience

Work experience is one of the most important variables in Singer and Willett
(2003) as it indicates time and makes longitudinal analysis possible. It is
desirable to calculate this rather than using the survey year for time because
it more accurately reflects a person’s time in the workforce. Thus, in the
spirit of refreshing the data to the newest round of the survey, this variable
needs to be calculated from other variables provided. It is not
straightforward. We start with the definition of experience in Singer and
Willett (2003).

Experience is years after entering the labor force. It represents the
difference between the day an individual enters the labor force
(EMPLOYERS_ALL_STARTDATE_ORIGINAL.01~Y_XRND in the database) relative to the
date of the survey, which we call yr_wforce. However, using this calculation
produced numbers that do not quite match the original data.

Reading the section titled “Topical Guide to the Data” in the guide suggests
it should be calculated based on the variable “number of weeks worked since
the last interview.” This would remove periods of unemployment which makes
sense when measuring experience while actually working. Since it is only
measured since the last interview, this needs to be cumulated for each survey
year. This produced results more similar to the original data (as discussed
further in Section 4.3).

###  3.3 Initial data analysis

According to Huebner, Vach, and Cessie (2016), initial data analysis (IDA) is
the step of inspecting and screening the data after collection to ensure the
data is clean, valid, and ready to be deployed in the later analyses. This is
supported by Chatfield (1985), who argues that the two main objectives of IDA
are data description, which is to assess the structure and the quality of the
data, and model formulation without any formal statistical inference.

In this paper, we conduct an IDA or a preliminary data analysis to assess the
validity of the variable values in the cohort of data that the NLSY provides.
The first step is validating numerical summaries of the raw data are the same
as reported by NLSY79. This is followed by graphical summaries using methods
available in ggplot2 (Wickham 2016) and brolgar (Tierney, Cook, and Prvan
2020).

The respondents’ ages ranged from 12 to 22 when first interviewed in 1979.
Hence, we validate whether all respondents were in this range in the data we
extracted. Additionally, the NLSY also provides the number of the survey
cohort by their sex (6,403 males and 6,283 females) and race (7,510 Non-
Black/Non-Hispanic; 3,174 Black; 2,002 Hispanic). To validate this, we used
the demog_nlsy79, i.e., the data with the survey years 1979 sample. Tables 1
and 2 suggest the demographic data we had is consistent with the sample
information in the database.

Figure 2: Longitudinal profiles of wages for a random sample of 36 individuals
in the pre-cleaned data. There is considerable variation in wages. Some
individuals (2799, 11041, 11146) are only measured for a short period. Some
individuals (8296, 9962) possibly have errors in wages in some years, because
of the extreme fluctuation.

In the next step, we explore the mean hourly wage data of samples of
individuals. The purpose is to examine the common patterns and check the
quality. A random sample of 36 individuals is chosen (using the sample_n_keys
function in brolgar). Their longitudinal profiles are plotted, faceted by id,
and using free $y$ scales so that the individual patterns can be examined
(Figure 2). There is a lot of variability from one individual to another and
substantial fluctuation in wages at different times for most individuals. Some
individuals (2799, 11041, 11146) are only measured for a short period. Some
individuals (8296, 9962) possibly have errors in wages in some years because
of the extreme fluctuation. These need to be inspected more closely. It is
also important to note that some shorter profiles indicate some individuals
have left the study before it has finished. Checking whether the demographics
of the early departing are similar to those who remain in the study is an
important part of any downstream analysis to account for the bias induced by
the inadvertent censoring.

Figure 3 shows an alternative way to check the data quality. Plot (A) is the
spaghetti plot where all profiles are shown, and it can be seen that there are
unbelievably high wage values (up to $60,000/hour) for some individuals,
mostly around 1990. Plot (B) shows side-by-side boxplots of the three number
summaries (minimum, median, and maximum) for all individuals. This tells us
there are a number of individuals with unbelievably high maximum wages. Plot
(C) shows the profile for an individual with a maximum wage that is not so
extreme but still indicates a problem: their wages are consistently low except
for one year where they earned close to $1200/hour. This does not seem to be
reasonable and leads us to use a procedure to detect and fix these temporal
anomalies.

Table 1: Frequency table of the age at the start of the survey in NSLY79 cohort in the extracted data. Age | Number of individuals
---|---
15 | 1,265
16 | 1,550
17 | 1,600
18 | 1,530
19 | 1,662
20 | 1,722
21 | 1,677
22 | 1,680
Table 2: Contingency table for sex and race for the extracted NLSY79 demographic data. The percentage (rounded to closest 1%) is out of the total corresponding to row. | Race |
---|---|---
Sex | Hispanic | Black | Non-Black, Non-Hispanic | Total
F | 1,561 (25%) | 1,002 (16%) | 3,720 (59%) | 6,283
M | 1,613 (25%) | 1,000 (16%) | 3,790 (59%) | 6,403
Total | 3,174 (25%) | 2,002 (16%) | 7,510 (59%) | 12,686
Figure 3: Summary plots to check the data after the tidying stage: (A)
longitudinal profiles of wages for all individuals 1979-2018, (B) boxplots of
minimum, median, and maximum wages of each individual, (C) and one individual
(id=39) with an unusual wage relative to their years of data. It reveals that
some values of hourly wages are unbelievable, and some individuals have
extremely unusual wages in some years. Accordingly, more cleaning is necessary
to treat these extreme values.

Extremely high values were also found in the total hours of work, where some
observations reported having worked for 420 hours a week in total. According
to Pergamit et al. (2001), one of the flaws of the NLSY79 employment data is
the NLSY79 collects the information of the working hours since the last
interview. Thus, it might be challenging for the respondents to track the
within-job hours’ changes between survey years, especially for the respondents
with fluctuating working hours or seasonal jobs. It even has been more
challenging since 1994, after which respondents were only surveyed every other
year and thus had to recall two full years’ job history. This shortcoming
might also contribute to the fluctuation of one’s wages data.

###  3.4 Replacing extreme values

A robust linear regression model using the rlm function from MASS package
(Venables and Ripley 2002) is used to treat the extreme values in the data.
The robustness weight is used to determine if a value should be replaced with
the fitted value from the model. This is constructed for each ID utilizing the
nest and map function from tidyr (Wickham 2020) and purrr (Henry and Wickham
2020), respectively. An alternative approach would be a robust linear mixed
model using robustlmm (Koller 2016).

Figure 4 shows the profiles for a sample of 30 individuals the mean hourly
wage before and after extreme values are replaced. The plot shows fluctuations
in wages remain, but the large spikes (in this sample, individuals 8296,
9962), which are considered implausible, are replaced.

The challenging part of detecting an anomaly using the robustness weight is
determining the weight threshold beyond which the observations are considered
outliers. It is also important not to be too aggressive in outlier removal,
because this would overly smooth the data, and ignore actual workforce
experiences important wage volatility such as, abnormally high wages for a
fleeting amount of time. To explore the risk of being overly vigorous in
labeling observations as outliers, a range of threshold value were assessed.
Numerous samples of individuals were drawn for each threshold, and the
smoothness of the profiles was examined. Overly smooth profiles would indicate
the replacement of values was too severe, resulting in removing the very
interesting volatility of wages seen in many individuals. This was achieved
using a purpose written shiny (Chang et al. 2020) app (provided with the code
of this paper). Observations where a change was made are tracked with a new
variable called is_pred, so that this effect in the downstream analyses can be
monitored.

Figure 4: Comparison between the original (black dots) and the corrected
(solid grey) mean hourly wage for same sample of individuals as shown in
Figure 2. A robust linear model prediction was used to identify and correct
mean hourly wages value. The extreme spikes, corresponding to implausible
wages, have been replaced with values more similar to wages in neighboring
years for individuals 8296 and 9962, but otherwise the profiles have not
changed.

Figure 5 shows the summary statistics after removing extremes. The highest
wage overall is now around $1000. Plot (A) shows a more reasonable spaghetti
plot, where there are some profiles with high wages, but most profiles have
wages under $300, and there has been a steady increase in wages over the
years. Plot (B) shows there are still a small number of individuals with high
maximum wages. Plot (C) shows the profile for ID=39 after imputing the extreme
value. The wages for this individual increase over the years, and do fluctuate
some between 1900 and 2005.

Figure 5: Re-make of the summary plots of the fully processed data suggest it
is now in a reasonable state: (A) longitudinal profiles of wages for all
individuals 1979-2018, (B) boxplots of minimum, median, (C) and maximum wages
of each individual, and one individual with an unusual wage relative to their
years of data.

###  3.5 Recap

There are many steps and decisions made to go from raw to input to valid data.
Figure 6 summarizes these in order to create a refreshed wages data set.

Figure 6: The stages of data cleaning from the raw data to get three datasets
contained in yowie. “# of individuals” means the number of respondents
included in each stage, while “# of observations” means the number of rows in
the data. The color represents the stage of data cleaning in the statistical
value chain (M. P. J. van der Loo and de Jonge 2021). Pink, blue, and green
represent the raw, input, and valid data, respectively.

The list of variables provided in the three new datasets are as follows:

demog_nlsy79 :

  1. 1.

id: A unique individual’s ID number.

  2. 2.

age_1979: The age of the individual in 1979.

  3. 3.

sex: Sex of the individual, f = Female and m = Male.

  4. 4.

race : Race of the individual belong to, NBH = Non-Black, Non-Hispanic; H =
Hispanic; B = Black.

  5. 5.

hgc: The highest grade completed ever.

  6. 6.

hgc_i: Integer value of the highest grade completed ever.

  7. 7.

hgc_1979: The highest grade completed in 1979 (integer value).

  8. 8.

ged: Whether the individual had a high school diploma or Graduate Equivalency
Degree (GED). 1: High school diploma; 2: GED; 3: Both.

wages and its subset for high school dropouts cohort wages_hs_do:

  1. 1.

id: A unique individual’s ID number. This is the key of the data as we saved
the data as a tsibble object.

  2. 2.

year: The year the observation was taken. This is the index of the data.

  3. 3.

wage: The mean of the hourly wages the individual gets at each of their
different jobs. The value could be a weighted or an arithmetic mean. The
weighted mean is used when the information of hours of work as the weight is
available. The mean hourly wage could also be a predicted value if the
original value is considered influential by the robust linear regression as
part of data cleaning.

  4. 4.

age_1979: The age of the individual in 1979.

  5. 5.

sex: Sex of the individual, f = Female and m = Male.

  6. 6.

race: Race of the individual belong to, NBH = Non-Black, Non-Hispanic; H =
Hispanic; B = Black.

  7. 7.

grade: Integer value of the highest grade completed corresponding to the
survey year.

  8. 8.

hgc: The highest grade completed ever.

  9. 9.

hgc_i: Integer value of the highest grade completed ever.

  10. 10.

hgc_1979: The highest grade completed in 1979 (integer value).

  11. 11.

ged: Whether the individual had a high school diploma or Graduate Equivalency
Degree (GED). 1: High school diploma; 2: GED; 3: Both.

  12. 12.

njobs: Number of jobs that an individual has.

  13. 13.

hours: The total number of hours the individual usually works per week.

  14. 14.

stwork: The year when the individual started to work.

  15. 15.

yr_wforce: The length of time in the workforce in years (year \- stwork).

  16. 16.

exp: Work experience, i.e., the number of years of working.

  17. 17.

is_wm: Whether the mean hourly wage is weighted mean, using the hour work as
the weight, or regular/arithmetic mean. TRUE = is weighted mean. FALSE = is
regular mean.

  18. 18.

is_pred: Whether the mean hourly wage is a predicted value of RLM or not.

##  4 Comparison of refreshed with the original data

The original set, containing wages of high school dropouts (Singer and Willett
2003) from 1979 through to 1994, is available in the R package brolgar
(Tierney, Cook, and Prvan 2020). To compare the refreshed data with the
original, a subset needs to be matched. There are numerous ways to do this,
with the simplest being to extract the individuals based on their id being
part of the original data and restricting the longitudinal measurements to the
same years. However, we decided to try to replicate the process, as suggested
by the description of the original data. This requires first identifying
individuals who dropped out of high school.

###  4.1 Filtering: Determining who is a dropout

There is no explicit explanation of how the dropouts cohort is determined in
the original data. Hence, we use the high school dropouts criteria from Wolpin
(2005), which are:

  1. 1.

An individual whose highest grade completed (hgc) is reported to be less than
12th grade, or

  2. 2.

An individual whose highest grade completed (hgc) is reported to be at least
12th grade and has received a GED (i.e. ged code is 2).

An additional criterion from Singer and Willett (2003) is to only include
males aged between 14 and 17 years old in 1979. With this filtering, we
obtained 670 individuals in the refreshed data compared to 888 individuals in
the original data. To investigate the reason for the difference, individuals
from the original and refreshed dataset were matched by id. This revealed
several reasons for the disparity:

  1. 1.

173 individuals were more than 17 years old in 1979. Thus, it looks like the
description of the original data is not quite accurate that there are people
older than 17 in the subset. Our decision is to also include them in the
refreshed data as the new data contains an age variable, so analysts could
filter them later.

  2. 2.

79 individuals were less than or equal to 17 years old in 1979. However, they
were not captured in the refreshed data because:

    1. i.

35 of them completed at least 12th grade with a diploma instead of GED (ged
variable is coded as 1). This suggests they are not dropouts, and so we
excluded them from the refreshed data.

    2. ii.

The information about ged is missing in 38 individuals. We decided to include
them in the refreshed data.

    3. iii.

3 individuals have both diploma and GED (ged is coded as 3). These were kept
in the refreshed data.

    4. iv.

12 individuals do not exist in wages data because they have participated in
less than 3 rounds of the survey.

The filtering was re-applied using these decisions, resulting in a refreshed
dropout subset containing 863 individuals.

###  4.2 Summaries of original with refreshed dropouts data

Because the original data does not have the year of collection, it is not
possible to merge the two subsets directly. Merging longitudinal data requires
both the key (id) and the index (ideally survey year). In the original data,
the experience variable is the time index, and it was not possible to exactly
match this for the refreshed data. Thus, comparisons of the two sets have to
be conducted in a two-sample fashion rather than a matched sample.

Figure 7 contains summaries of corresponding variables in the two subsets.
Plot (A) shows a back-to-back bar chart of the highest grade completed. The
two sets are almost the same, but small differences remain. Plots (B) and (C)
show stacked density plots of experience and log wages, respectively. The
distributions are relatively close, with more differences in wages as would be
expected because the refreshed data is not inflation-adjusted.

Figure 7: Comparison of original and refreshed data: (A) highest grade
completed, (B) experience and (C) log wages. Some difference in wages would be
expected because the refreshed data is not inflation-adjusted, but the two
sets are reasonably similar.

###  4.3 The takeaways

There are several aspects of the original data that were difficult to
replicate. The calculation of the work experience is not clearly articulated.
Singer and Willett (2003) describe the temporal variable as the years of
experience since the first day of work. This variable is not explicitly
available in the database. From the NLSY79 topical guide (Bureau of Labor
Statistics, U.S. Department of Labor 2021b), we find several variables are
tagged as work experience-related variables. One of them is the weeks worked
since the last interview. This is used to calculate the variable. It produces
reasonably similar but not exactly the same values. Because the original data
set did not include the year of the survey, it cannot be precisely compared.

The highest grade completed has some confusion. There are several ways this is
reported, including the highest grade ever completed and also the hgc at each
survey year. To match the original data, it is appropriate to use the hgc
while in high school. The documentation suggests this variable is available,
but it is not actually present in the data. Hence, to match the original data,
we have calculated the hgc to match the hgc achieved in the years between 1979
and 1994 based on the yearly survey value. The result does not exactly match
the original for a few individuals.

In the original dataset, wages were inflation-adjusted to 1990 prices. This is
not done for the refreshed data because we plan to keep refreshing it as the
data is added to and released from the survey. Instead, we have provided a
function in the R package, yowie, for users to conduct the inflation
adjustment when they are ready to analyze the data.

It is important to note that the treatment of unlikely wages differs in the
refreshed data. In the original data by Singer and Willett (2003), wages
greater than $75 are set to be missing. However, this value is too low to be
set as the maximum threshold, and it doesn’t take into account temporal
neighbors for an individual. We opted to use the weights from a robust linear
regression to determine what should be regarded as extreme and imputed them
with their predicted values as described in Section 3.4.

Lastly, the matching of individuals in the dropouts subset of the refreshed
data with the original data was done using their ids. (Only males were in the
original data.) It was comforting to see that all of the individuals in both
sets do match on sex (all are male!), and race.

##  5 Summary

Any longitudinal dataset used for education should have a sufficiently
reproducible process to be refreshed with new data. The NLSY79 dataset is a
great teaching tool that has become outdated, both because the dataset stops
in 1994, and because the demographic data could be handled more delicately.
This paper has illustrated the necessary steps and decisions made to take a
particular open data set and make it a textbook data set, ready for the
classroom or research. In the first stage, we showed the steps performed to
get the data from the NLSY79 database. The data format was converted to a tidy
format for more flexibility in cleaning and exploring. Initial data analysis
was conducted to investigate and screen the quality of the data. We found and
provided a fix for many anomalous observations in wages using a robust linear
regression model. The refreshed data is compared with the original set using a
variety of numerical summaries and graphics. The current subset is made
available in a new R package, called yowie.

The data cleaning process is documented, and the code has been made available.
These provide the opportunity to again refresh the textbook data as new data
is published into the NLSY79 database. Determining an appropriate robustness
weight from which to threshold unusual observations was conducted using a
shiny (Chang et al. 2020) app, which is provided with the code, and the
choices used in the refreshed data are documented.

Various difficulties were encountered in trying to refresh the data, which
include:

  * •

Determining which records should be downloaded from the database.

  * •

Calculating experience in the workforce requires comparing the date of the
first job with the first year the individual was recorded, both of which are
available in the database.

  * •

Treating the extreme values since there are many unusually high hourly wages,
e.g., greater than $60,000 per hour.

  * •

Determining the dropouts subset as there is no explicit variable in the
database recording high school dropout, which means we needed to compare the
date of 12th grade with their GED status.

  * •

Matching IDs from the original data with those in the refreshed data do refer
to the same person based on the demographic information available.

Ultimately, the refreshed data is reasonably similar to the original but
unsatisfactorily far from it. The last step required would be to inflation-
adjust wages. This is better to do with each wave of new data added so that it
is relative to the last date in the data. Our decision was to provide the raw
wages and include code to make the adjustment as part of the package.

Some readers may disagree with our decisions made to produce the refreshed
textbook data and may have better insight than us in producing more
appropriate textbook data. We do not assert that we have produced the best
textbook data, but rather we describe our journey to provide a reasonable
textbook data set. All code and documentation are provided for transparency.
Readers could use this to make different decisions or provide suggestions
through the package for better choices. Future updates of the yowie package
may contain additional variables, or filters of the full set if it is deemed
important.

For the data providers, we recommend a better validation system with clear
rules applied at data entry and alternative output formats, such as a tidy
format, which would help users better use their resources. The problem with
many wages records is that there are implausible values or confusion on how to
record wages for multiple jobs. These values can be validated with simple
checks at data entry. Providing an open data resource is also accompanied by
the responsibility that the data, especially as valuable as this, is reliable.
Users need to be able to trust the data.

Why is the work presented here useful for teachers and students? There are
several reasons:

  * •

It primarily illustrates the steps in cleaning and processing a well-known
longitudinal data so that it can be further refreshed as a textbook data set
as new data emerges. Choices made during the cleaning and processing can
affect findings made with the data, are transparent, and can be a basis for
classroom discussions.

  * •

Some of the methods shown, and lessons learned about the data, can also be
considered generally for data cleaning, initial data analysis and exploratory
data analysis, useful for working with other data sets.

  * •

The refreshed data can be used for teaching advanced linear modeling, using
almost identical code from Singer and Willett (2003), and thus make those
lessons more current for today’s students.

  * •

It provides the data to develop a case study for teaching exploratory data
analysis of longitudinal data, focusing on the individual experience. The
vignettes of the brolgar (Tierney, Cook, and Prvan 2020) package provide
useful guides for the case study content.

The wages data provides a good opportunity to discuss the difference between
statistics for public policy and statistics that relate to the individual.
Public policy is based on models, yielding averages that might vary across
strata. Modeling the wages relative to workforce experience with demographic
covariates, we learn that there are significantly different patterns. More
education leads to increasingly higher wages, which is a satisfying result for
educators. It means education makes a difference in the wage experience, so
public policy encouraging education is a data-supported action. We would also
learn, although not presented here, that wages on average differ by race
category. This is a disturbing finding because there is no rationale for such
a difference in a fair society. This would provide support for action in
public policy to remove this overall average effect. It also provides an
example for educators to explain the use of data to support public policy
action.

On an individual level, people might want to see how their characteristics
relate to those in the dataset. To do this, the individual profiles need to be
explored. Predominantly, we would learn the variation from one individual to
another is far more than the variation between demographic strata. For
example, many individuals with lower educational attainment earn very high
wages. From a statistics and data science educator perspective, more focus and
more methodology for this type of statistics need to be included in the
curriculum.

The above interpretations of results from analyzing the wages data rely on
trusting that the data provided is accurate and valid. The wages data is
collected by a reputable organization, but we found the data has obvious
errors that should be corrected. This paper has illustrated procedures and
guidelines to achieve valid data and provides the code and details for it to
be reproduced and modified if deemed appropriate.

Having trustworthy data is imperative for statistics and data science
education. Whenever one uses a textbook data set that is perceived as relating
to the students’ lives, there will be interpretations made. The wages data is
an example of this. Students can take away multiple interpretations from the
data (e.g. wages increase with experience, education, and race) based on the
teaching focus to learn about the world. If one uses data examples that are
synthetic, instilled with our own inherent prejudices (e.g., sex and race), or
data that has been poorly processed containing errors, we, as educators, are
being irresponsible because the societal message taught to students may be
flawed. This paper demonstrates the process of producing a trustworthy data
set for teaching.

##  6 Acknowledgements

We would like to thank Aarathy Babu for the insight and discussion during the
writing of this paper.

The entire analysis is conducted using R (R Core Team 2020) in RStudio IDE
using these packages: tidyverse (Wickham et al. 2019), ggplot2 (Wickham 2016),
dplyr (Wickham et al. 2020), readr (Wickham and Hester 2020), tidyr (Wickham
2020), stringr (Wickham 2019), purrr (Henry and Wickham 2020), brolgar
(Tierney, Cook, and Prvan 2020), patchwork (Pedersen 2020), kableExtra (Zhu
2019), MASS (Venables and Ripley 2002), janitor (Firke 2020), and tsibble
(Wang, Cook, and Hyndman 2020). The paper was generated using knitr (Xie 2014)
and rmarkdown (Xie, Dervieux, and Riederer 2020).

##  7 Supplementary Materials

  * •

Code: R script to reproduce data tidying and cleaning is available at
https://numbats.github.io/yowie/articles/process-data.html. The code for the
extreme value handling with robust linear fits is in
https://numbats.github.io/yowie/articles/input-to-valid-data.html.

  * •

R Package: yowie is an R data package that contains 3 datasets, namely the
high school mean hourly wage data, high school dropouts mean hourly wage data,
and demographic data of the NLSY79 cohort. This package can be accessed from
https://github.com/numbats/yowie.

  * •

shiny app: An interactive shiny web app to visualise the effect of selecting
different weight threshold for substituting the wages data to its predicted
value from a fit of the robust linear regression model. This app can be
accessed at https://ebsmonash.shinyapps.io/yowie_app/ with the source code
provided https://github.com/numbats/yowie/tree/master/app.

##  8 Data Availability Statement

The authors confirm that the data supporting the findings of this study are
available within the supplementary materials.

## References

reAnderson, Edgar. 1935. “The Irises of the Gaspe Peninsula.” _Bulletin of the
American Iris Society_ 59: 2–5.

preBureau of Labor Statistics, U.S. Department of Labor. 2021a. “National
Longitudinal Survey of Youth 1979 Cohort, 1979-2016 (Rounds 1-28).” Produced
and distributed by the Center for Human Resource Research (CHRR), The Ohio
State University. Columbus, OH, through https://www.nlsinfo.org/bibliography-
citing-nls-data.

pre———. 2021b. “National Longitudinal Survey of Youth 1979 Cohort, Topical
Guide to the Data.” https://www.nlsinfo.org/content/cohorts/nlsy79/topical-
guide/employment/work-experience.

preChang, Winston, Joe Cheng, JJ Allaire, Yihui Xie, and Jonathan McPherson.
2020. _shiny: Web Application Framework for R_.
https://CRAN.R-project.org/package=shiny.

preChatfield, C. 1985. “The Initial Examination of Data.” _Journal of the
Royal Statistical Society. Series A. General_ 148 (3): 214–53.

preCooksey, Elizabeth C. 2017. “Using the National Longitudinal Surveys of
Youth (NLSY) to Conduct Life Course Analyses.” In _Handbook of Life Course
Health Development_ , edited by Richard M. Lerner Neal Halfon Christoper B.
Forrest, 561–77. Cham: Springer.
https://doi.org/https://doi.org/10.1007/978-3-319-47143-3_23.

preDasu, Tamraparni, and Theodore Johnson. 2003. _Exploratory Data Mining and
Data Cleaning_. Wiley Series in Probability and Statistics. Hoboken: WILEY.

preFirke, Sam. 2020. _janitor: Simple Tools for Examining and Cleaning Dirty
Data_. https://CRAN.R-project.org/package=janitor.

preFullilove, M. T. 1998. “Comment: Abandoning “Race" as a Variable in Public
Health Research–an Idea Whose Time Has Come.” _American Journal of Public
Health_ 88 (9): 1297–98.

preGrimshaw, Scott D. 2015. “A Framework for Infusing Authentic Data
Experiences Within Statistics Courses.” _The American Statistician_ 69 (4):
307–14. https://doi.org/https://doi.org/10.1080/00031305.2015.1081106.

preHeidari, S., T. F. Babor, P. De Castro, S. Tort, and M. Curno. 2016. “Sex
and Gender Equity in Research: Rationale for the SAGER Guidelines and
Recommended Use.” _Research Integrity and Peer Review_ 1 (2).
https://doi.org/https://doi.org/10.1186/s41073-016-0007-6.

preHenry, Lionel, and Hadley Wickham. 2020. _purrr: Functional Programming
Tools_. https://CRAN.R-project.org/package=purrr.

preHorst, Allison Marie, Alison Presmanes Hill, and Kristen B Gorman. 2020.
_Palmerpenguins: Palmer Archipelago (Antarctica) Penguin Data_.
https://doi.org/10.5281/zenodo.3960218.

preHuebner, Marianne, Werner Vach, and Saskia le Cessie. 2016. “A Systematic
Approach to Initial Data Analysis Is Good Research Practice.” _The Journal of
Thoracic and Cardiovascular Surgery_ 151 (1): 25–27.

preIlk, Ozlem. 2004. “Exploratory Multivariate Longitudinal Data Analysis and
Models for Multivariate Longitudinal Binary Data.” PhD thesis, Iowa State
University. https://doi.org/10.31274/rtd-180813-11012.

preKennedy, Lauren, Katharine Khanna, Daniel Simpson, and Andrew Gelman. 2020.
“Using Sex and Gender in Survey Adjustment.” https://arxiv.org/abs/2009.14401.

preKim, A. Y, C. Ismay, and J. Chunn. 2018. “The Fivethirtyeight R Package:
“Tame Data" Principles for Introductory Statistics and Data Science Courses.”
_Technology Innovations in Statistics Education_ 11 (1).
https://doi.org/https://doi.org/10.5070/T511103589.

preKoller, Manuel. 2016. “robustlmm: An R Package for Robust Estimation of
Linear Mixed-Effects Models.” _Journal of Statistical Software_ 75 (6): 1–24.

preMoncrief, Marc. 2015. “By the Numbers - the Average Australian Doesn’t
Exist … Not a Single One of Us Is ’Normal’.” https://bit.ly/smh-not-normal.

preOffice of Management and Budget. 1997. “Revisions to the Standards for the
Classification of Federal Data on Race and Ethnicity.”
https://www.govinfo.gov/content/pkg/FR-1997-10-30/pdf/97-28653.pdf.

preOpen Knowledge Foundation. 2021. “Open Definition. Defining Open in Open
Data, Open Content, and Open Knowledge.” 2021.
http://opendefinition.org/od/2.1/en/.

prePedersen, Thomas Lin. 2020. _patchwork: The Composer of Plots_.
https://CRAN.R-project.org/package=patchwork.

prePergamit, Michael R., Charles R. Pierret, Donna S. Rothstein, and Jonathan
R. Veum. 2001. “Data Watch: The National Longitudinal Surveys.” _The Journal
of Economic Perspectives_ 15 (2): 239–53.

preR Core Team. 2020. _R: A Language and Environment for Statistical
Computing_. Vienna, Austria: R Foundation for Statistical Computing.
https://www.R-project.org/.

preSinger, Judith D, and John B Willett. 2003. _Applied Longitudinal Data
Analysis: Modeling Change and Event Occurrence_. Oxford u.a: Oxford Univ. Pr.

preStodel, Megan. 2020. “Stop Using Iris.”
https://www.meganstodel.com/posts/no-to-iris/.

preTierney, Nicholas, Di Cook, and Tania Prvan. 2020. _brolgar: BRowse Over
Longitudinal data Graphically and Analytically in R_.
https://github.com/njtierney/brolgar.

preTukey, John W. (John Wilder). 1977. _Exploratory Data Analysis_. Addison-
Wesley Series in Behavioral Science. Reading, Mass.: Addison-Wesley Pub. Co.

prevan der Loo, Mark P. J., and Edwin de Jonge. 2021. “Data Validation
Infrastructure for R.” _Journal of Statistical Software_ 97 (10): 1–31.
https://doi.org/10.18637/jss.v097.i10.

prevan der Loo, Mark, and Edwin de Jonge. 2018. _Statistical Data Cleaning
with Applications in r_.

preVenables, W. N., and B. D. Ripley. 2002. _Modern Applied Statistics with
s_. Fourth. New York: Springer. http://www.stats.ox.ac.uk/pub/MASS4.

preWang, Earo, Dianne Cook, and Rob J Hyndman. 2020. “A New Tidy Data
Structure to Support Exploration and Modeling of Temporal Data.” _Journal of
Computational and Graphical Statistics_ 29 (3): 466–78.
https://doi.org/10.1080/10618600.2019.1695624.

preWickham, Hadley. 2011. “The Split-Apply-Combine Strategy for Data
Analysis.” _Journal of Statistical Software, Articles_ 40 (1): 1–29.
https://doi.org/https://doi.org/10.18637/jss.v040.i01.

pre———. 2014. “Tidy Data.” _Journal of Statistical Software_ 59 (10): 1–23.

pre———. 2016. _ggplot2: Elegant Graphics for Data Analysis_. Springer-Verlag
New York. https://ggplot2.tidyverse.org.

pre———. 2019. _stringr: Simple, Consistent Wrappers for Common String
Operations_. https://CRAN.R-project.org/package=stringr.

pre———. 2020. _tidyr: Tidy Messy Data_.
https://CRAN.R-project.org/package=tidyr.

preWickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy
D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome
to the tidyverse.” _Journal of Open Source Software_ 4 (43): 1686.
https://doi.org/10.21105/joss.01686.

preWickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2020.
_dplyr: A Grammar of Data Manipulation_.
https://CRAN.R-project.org/package=dplyr.

preWickham, Hadley, and Jim Hester. 2020. _readr: Read Rectangular Text Data_.
https://CRAN.R-project.org/package=readr.

preWolpin, Kenneth I. 2005. “National Longitudinal Survey of Youth 1979
Cohort, 1979-2016 (Rounds 1-28).” Published by Bureau of Labor Statistics,
U.S. Department of Labor. https://www.bls.gov/opub/mlr/2005/02/art3full.pdf.

preXie, Yihui. 2014. “Knitr: A Comprehensive Tool for Reproducible Research in
R.” In _Implementing Reproducible Computational Research_ , edited by Victoria
Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.
http://www.crcpress.com/product/isbn/9781466561595.

preXie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. _R Markdown
Cookbook_. Boca Raton, Florida: Chapman; Hall/CRC.
https://bookdown.org/yihui/rmarkdown-cookbook.

preZhu, Hao. 2019. _kableExtra: Construct Complex Table with ’kable’ and Pipe
Syntax_. https://CRAN.R-project.org/package=kableExtra.

p


# Data-Driven Inverse of Linear Systems
and Application to Disturbance Observers

Yongsoon Eun, Jaeho Lee, and Hyungbo Shim This work was supported by the
National Research Foundation of Korea(NRF) grant funded by the Korea
government(MSIT)(NRF2022R1A2C1010777).Yongsoon Eun and Jaeho Lee are with the
Department of Electrical Engineering and Computer Science, DGIST, Daegu 42988,
South Korea. {yeun, jaeho.lee} @ dgist.ac.krHyungbo Shim is with ASRI,
Department of Electrical and Computer Engineering, Seoul National University,
Seoul, Korea<EMAIL_ADDRESS>
###### Abstract

This work develops a data-based construction of inverse dynamics for LTI
systems. Specifically, the problem addressed here is to find an input sequence
from the corresponding output sequence based on pre-collected input and output
data. The problem can be considered as a reverse of the recent use of
behavioral approach, in which the output sequence is obtained for a given
input sequence by solving an equation formed by pre-collected data. The
condition under which the problem gives a solution is investigated and turns
out to be $L$-delay invertibility of the plant and a certain degree of
persistent excitation of the data input. The result is applied to form a data-
driven disturbance observer. The plant dynamics augmented by the data-driven
disturbance observer exhibits disturbance rejection without the model
knowledge of the plant.

##  I Introduction

Seminal work of the behavioral approach [1] for dynamic systems has been
receiving growing attention in recent years [2, 3, 4]. Instead of using
matrices $A$, $B$, $C$ and $D$ for linear dynamics, which are often obtained
from the first principle on the target dynamics, the behavioral approach uses
a set of equations to determine if a pair of input and output belongs to the
target dynamics. The equations involve Hankel matrices for linear dynamics,
which are obtained by input/output trajectories that satisfy a condition of
persistently exciting. This approach suits when data are available but
applying the first principle for modeling is not straightforward. Recent
advances in this area include [5, 6, 3, 7].

Particularly, as given in [8], the fundamental result of this approach is that
we can determine a future output for a given input using Hankel matrices built
by using previously collected input and output data. No identification of
system matrices (typically $A$, $B$, $C$, $D$) is necessary. A question that
this paper is concerned is the inverse of what is just described: can we
determine the input that generated a given output using previously collected
input and output data? In other words, the question pertains to finding
conditions and methods how we can build an inverse dynamics of the given
system using collected data. If this is successfully carried out, an immediate
application is to build a disturbance observer [9, 10] with data, that
identifies the disturbance that affected the output. Another application is to
find a suitable constrained input that yields an output close to the desired
output under constraints.

Similar question has been asked and answered in [11, 8]. However, the result
is restricted to the case where $D$ is invertible. This work does not assume
the invertibility of $D$. Instead the notion of $L$-delay invertibility is
invoked from the literature [12].

This paper is organized as follows: Section II provides a brief review of the
invertibility of discrete-time LTI systems. Section III gives main results. An
application of disturbance observer is given in Section IV. Section V
concludes the paper.

##  II Review of LTI System Invertibility

Consider a discrete-time linear-time-invariant system

$\displaystyle x(k+1)$ $\displaystyle=Ax(k)+Bu(k)$ (1a) $\displaystyle y(k)$
$\displaystyle=Cx(k)+Du(k)$ (1b)

where $x\in\mathbb{R}^{n}$ is the state, $u\in\mathbb{R}^{m}$ is the input,
$y\in\mathbb{R}^{p}$ is the output, and $k\in\mathbb{Z}$ is the discrete time
index. Then the transfer function matrix of given system (1) is given by
$G(z)=C(zI_{n}-A)^{-1}B+D$.

######  Definition 1 (see [12]).

A proper transfer function matrix $\hat{G}(z)$ is an $L$-delay inverse system
for given system (1), where $L$ is a nonnegative integer, if

$\displaystyle\hat{G}(z)G(z)=\dfrac{1}{z^{L}}I_{m}.$ (2)

Note that $G(z)$ is an $p\times m$ transfer function matrix and $\hat{G}(z)$
is a $m\times p$ transfer function matrix. Thus, a necessary condition for (2)
is that $m\leq p$.

######  Definition 2 (see [12]).

The system (1) is invertible if it has an $L$-delay inverse for some finite
$L$. The least integer $L$ for which an $L$-delay inverse exists will be
called the inherent delay of the invertible system and is denoted by $L_{0}$.

For single input and single output systems of (1), the inherent delay $L_{0}$
is the relative degree.

Define, for an integer $t>0$,

$\displaystyle\begin{split}&\mathcal{T}_{t}\\\ &\mathrel{\mathop{\mathchar
58\relax}}=\begin{bmatrix}D&\mymathbb{0}&\mymathbb{0}&\cdots&\mymathbb{0}\\\
CB&D&\mymathbb{0}&\cdots&\mymathbb{0}\\\ CAB&CB&D&\cdots&\mymathbb{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
CA^{t-1}B&CA^{t-2}B&CA^{t-3}B&\cdots&D\end{bmatrix}.\end{split}$ (3)

In addition, define $\mathcal{T}_{0}=D$. Then, the following theorem gives the
necessary and sufficient condition for the existence of an $L$-delay inverse
for the system of (1).

######  Theorem 1 (see [12]).

The following are equivalent:

  1. 1.

The system (1) has an $L$-delay inverse.

  2. 2.

There exists $L$ such that

$\displaystyle{\rm rank}(\mathcal{T}_{L})-{\rm rank}(\mathcal{T}_{L-1})=m,$
(4)

where ${\rm rank}(\mathcal{T}_{-1})=0$.

  3. 3.

The rank of the matrix $G(z)$ is $m$.

As a result of Theorem 1, if the system has $L$-delay inverse, then the first
$m$ columns of $\mathcal{T}_{L}$ are independent.

##  III Data-based Construction of Inverse Systems

We introduce notations commonly used in the literature for data-driven system
and control. For a given signal $f(k)\in\mathbb{R}^{q}$, $k\in\mathbb{Z}$, and
$T\in\mathbb{N}$, let $f_{[k,k+T]}$ be defined by

$f_{[k,k+T]}\mathrel{\mathop{\mathchar 58\relax}}=\begin{bmatrix}f(k)\\\
f(k+1)\\\ \vdots\\\ f(k+T)\end{bmatrix}\quad\in\mathbb{R}^{q(T+1)}.$

Using the vectors $u_{[k,k+L]}$ and $y_{[k,k+L]}$, the system of (1) can be
written as

$\displaystyle y_{[k,k+L]}=\mathcal{O}_{L}x(k)+\mathcal{T}_{L}u_{[k,k+L]}$ (5)

where the matrix $\mathcal{O}_{L}$ is defined by

$\displaystyle\mathcal{O}_{L}=\begin{bmatrix}C^{T}\ (CA)^{T}\
\cdots(CA^{L})^{T}\end{bmatrix}^{T}.$ (6)

When the $L$-delay inverse exists for the system of (1), there exists an
$m\times p(L+1)$ matrix $\mathcal{K}$ satisfying

$\mathcal{K}\mathcal{T}_{L}=\begin{bmatrix}I_{m}&\mymathbb{0}_{m\times
mL}\end{bmatrix}.$

Such a matrix $\mathcal{K}$ exists due to Rouché-Capelli theorem [13] with the
fact that

${\rm rank}(\mathcal{T}_{L})={\rm rank}\left(\begin{bmatrix}\mathcal{T}_{L}\\\
I_{m}\ \mymathbb{0}_{m\times mL}\end{bmatrix}\right)$

which follows from the fact that the first $m$ columns of $\mathcal{T}_{L}$
are independent. Pre-multiplying $\mathcal{K}$ to both sides of (5) yields

$\displaystyle u(k)=-\mathcal{K}\mathcal{O}_{L}x(k)+\mathcal{K}y_{[k,k+L]},$
(7)

which represents the input $u(k)$ by the state and sequence of outputs. State
update equation with respect to $y_{[k,k+L]}$ is obtained by substituting (7)
in (1a). Then, the following dynamics is obtained:

$\displaystyle x(k+1)$ $\displaystyle=\tilde{A}x(k)+\tilde{B}y_{[k,k+L]}$ (8a)
$\displaystyle u(k)$ $\displaystyle=\tilde{C}x(k)+\tilde{D}y_{[k,k+L]}$ (8b)

where $\tilde{A}=A-B\mathcal{K}\mathcal{O}_{L}$, $\tilde{B}=B\mathcal{K}$,
$\tilde{C}=-\mathcal{K}\mathcal{O}_{L}$, and $\tilde{D}=\mathcal{K}$.

In order to provide more specific arguments, we first remind readers of Hankel
matrix and the notion of persistent excitation (PE) for a signal
$f(k)\in\mathbb{R}^{q}$. Denote by $\mathcal{H}_{t}(f_{[i,j]})$ the
$qt\times(j-i-t+2)$ dimensional Hankel matrix constructed using the values of
$f(i),\cdots,f(j)$:

$\displaystyle\mathcal{H}_{t}(f_{[i,j]})$
$\displaystyle=\begin{bmatrix}f(i)&f(i+1)&\cdots&f(j-t+1)\\\
f(i+1)&f(i+2)&\cdots&f(j-t+2)\\\ \vdots&\vdots&\ddots&\vdots\\\
f(i+t-1)&f(i+t)&\cdots&f(j)\\\ \end{bmatrix}.$

The notion of persistency of excitation of a segment from a signal is defined
as follows.

######  Definition 3 (see [1]).

The vector $f_{[0,T-1]}\in\mathbb{R}^{qT}$ from signal $f(k)\in\mathbb{R}^{q}$
is persistently exciting of order $t$ if the matrix
$\mathcal{H}_{t}(f_{[0,T-1]})\in\mathbb{R}^{qt\times(T-t+1)}$ has full row
rank.

As hinted by (8), data-based construction of inverse dynamics may involve
input and output sequences of different length.

######  Lemma 1.

A pair of vectors $u_{[0,T-1]}\in\mathbb{R}^{mT}$ and
$y_{[0,T-1+L]}\in\mathbb{R}^{p(T+L)}$ is a trajectory of the system (1) if and
only if there exists $u_{[T,T-1+L]}$ such that the pair $u_{[0,T-1+L]}$ and
$y_{[0,T-1+L]}$ is a length $T+L$ trajectory of the system (1).

Now, for some $T>0$ and $L\geq 0$, denote the input and output data collected
from an experiment for the system (1) by $u^{d}_{[0,T+L-1]}$ and
$y^{d}_{[0,T+L-1]}$, respectively.

######  Lemma 2.

Let the system (1) be controllable. Assume that $u^{d}_{[0,T-1]}$ is
persistently exciting of order $n+t+L$. Then, for any $t_{0}$, a pair of a
length $t$ input $u_{[t_{0},t_{0}+t-1]}$ and a length $t+L$ output
$y_{[t_{0},t_{0}+t+L-1]}$ is a trajectory of (1) if and only if there is
$g\in\mathbb{R}^{T-t+1}$ such that

$\displaystyle\begin{bmatrix}u_{[t_{0},t_{0}+t-1]}\\\
y_{[t_{0},t_{0}+t+L-1]}\end{bmatrix}=\begin{bmatrix}\mathcal{H}_{t}(u^{d}_{[0,T-1]})\\\
\mathcal{H}_{t+L}(y^{d}_{[0,T+L-1]})\end{bmatrix}g.$ (9)

###### Proof.

Let $u_{[t_{0},t_{0}+t-1]}$ and $y_{[t_{0},t_{0}+t+L-1]}$ be a trajectory of
(1). Then, by Lemma 1, there exists $u_{[t_{0}+t,t_{0}+t+L-1]}$ such that the
pair $u_{[t_{0},t_{0}+t+L-1]}$ and $y_{[t_{0},t_{0}+t+L-1]}$ is a trajectory
of (1). Then, by [5, Theorem 1] there exists $g$ such that

$\displaystyle\begin{bmatrix}u_{[t_{0},t_{0}+t+L-1]}\\\
y_{[t_{0},t_{0}+t+L-1]}\end{bmatrix}=\begin{bmatrix}\mathcal{H}_{t+L}(u^{d}_{[0,T+L-1]})\\\
\mathcal{H}_{t+L}(y^{d}_{[0,T+L-1]})\end{bmatrix}g.$ (10)

Hence (9) is satisfied. Conversely, let $u_{[t_{0},t_{0}+t-1]}$ and
$y_{[t_{0},t_{0}+t+L-1]}$ satisfy (9) for some $g$. Set the subsequent
sequence of $u$ as

$u_{[t_{0}+t,t_{0}+t+L-1]}=\mathcal{H}_{L}(u^{d}_{[t,T+L-1]})g.$ (11)

Then, $u_{[t_{0},t_{0}+t+L-1]}$ and $y_{[t_{0},t_{0}+t+L-1]}$ satisfy (10).
Again due to [5, Theorem 1] the pair is the trajectory of (1), and due to
Lemma 1, the pair $u_{[t_{0},t_{0}+t-1]}$ and $y_{[t_{0},t_{0}+t+L-1]}$ is a
trajectory of (1). ∎

Denote by $U_{p}$, $U_{f}$, $Y_{p}$ and $Y_{f_{L}}$ the following four Hankel
matrices formed by the data $u^{d}_{[0,T+L-1]}$, $y^{d}_{[0,T+L-1]}$ for some
$T_{p}>0$ and $T_{f}>0$ satisfying $T>T_{p}+T_{f}-1$

$\displaystyle U_{p}$
$\displaystyle=\mathcal{H}_{T_{p}}(u^{d}_{[0,T-T_{f}-1]}),\quad
U_{f}=\mathcal{H}_{T_{f}}(u^{d}_{[T_{p},T-1]}),$ $\displaystyle Y_{p}$
$\displaystyle=\mathcal{H}_{T_{p}}(y^{d}_{[0,T-T_{f}-1]}),\quad
Y_{f_{L}}=\mathcal{H}_{T_{f}+L}(y^{d}_{[T_{p},T+L-1]}).$

######  Lemma 3.

Let the system of (1) be observable and have an $L$-delay inverse with $L\geq
0$. Assume that $T_{p}$ be greater than or equal to the observability index of
(1). Then,

$\mathcal{N}\left(\begin{bmatrix}U_{p}\\\ Y_{p}\\\ U_{f}\\\
Y_{f_{L}}\end{bmatrix}\right)=\mathcal{N}\left(\begin{bmatrix}U_{p}\\\
Y_{p}\\\ Y_{f_{L}}\end{bmatrix}\right)$ (12)

where $\mathcal{N}(A)$ is the null space of matrix $A$.

###### Proof.

Denote $[x(T_{p}),x(T_{p}+1),\cdots,x(T-T_{f})]$ by $X$, and define
$v(k)=y^{d}_{[k,k+L]}$ for the simplicity of derivation. First, we show that
$U_{f}$ is determined by $X$ and $\mathcal{H}_{T_{f}}(v_{[T_{p},T-1]})$ by (7)
and (8) as

$\displaystyle U_{f}$ $\displaystyle=\begin{bmatrix}\tilde{C}\\\
\tilde{C}\tilde{A}\\\ \vdots\\\ \tilde{C}\tilde{A}^{T_{f}-1}\end{bmatrix}X$
$\displaystyle+\begin{bmatrix}\tilde{D}&\mymathbb{0}_{m\times
p(L+1)}&\cdots&\mymathbb{0}_{m\times p(L+1)}\\\
\tilde{C}\tilde{B}&\tilde{D}&\cdots&\mymathbb{0}_{m\times p(L+1)}\\\
\vdots&\ddots&\ddots&\vdots\\\
\tilde{C}\tilde{A}^{T_{f}-2}\tilde{B}&\tilde{C}\tilde{A}^{T_{f}-3}\tilde{B}&\cdots&\tilde{D}\end{bmatrix}$
$\displaystyle\quad\times\mathcal{H}_{T_{f}}(v_{[T_{p},T-1]}).$

Again using (3), (6), and (8), $X$ is determined by $U_{p}$ and $Y_{p}$ as

$\displaystyle
X=A^{T_{p}}\mathcal{O}_{T_{p}-1}^{\dagger}(Y_{p}-\mathcal{T}_{T_{p}-1}U_{p})+\mathcal{C}_{T_{p}-1}U_{p}$
(13)

where

$\displaystyle\mathcal{C}_{T_{p}-1}=\begin{bmatrix}A^{T_{p}-1}B&A^{T_{p}-2}B&\cdots&AB&B\end{bmatrix}$

and $\dagger$ implies the left-inverse of a matrix. Finally, let
$M_{i}\in\mathbb{R}^{p(L+1)\times p(T_{f}+L)}$ be a matrix that include the
identity matrix of size $p(L+1)$ from its $p(i-1)+1$-th column and all other
elements are zero. Then, it is easily seen that

$\displaystyle\mathcal{H}_{T_{f}}(v_{[T_{p},T-1]})=MY_{f_{L}}$ (14)

where $M^{T}=[M_{1}^{T},\cdots,M_{T_{f}}^{T}]$. Therefore,
$\mathcal{H}_{T_{f}}(v_{[T_{p},T-1]})$ is determined by $Y_{f_{L}}$. Combining
(13) and (14), we arrive that $U_{f}$ linearly depends on $U_{p}$, $Y_{p}$ and
$Y_{f_{L}}$. Therefore, (12) follows. ∎

Now we state the main result of this paper.

######  Theorem 2.

Let the system of (1) be controllable, observable and have an $L$-delay
inverse with $L\geq 0$. Assume that $T_{p}$ be greater than or equal to the
observability index of (1) and $u^{d}_{[0,T-1]}$ is PE of order
$n+T_{p}+T_{f}+L$. Then, for any trajectory $u_{[t_{0},t_{0}+T_{p}-1]}$ and
$y_{[t_{0},t_{0}+T_{p}+T_{f}+L-1]}$ of (1) where $t_{0}$ is arbitrary, the
input sequence $u_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}-1]}$ is uniquely determined
by

$u_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}-1]}=U_{f}g,$ (15)

where $g$ is a solution of

$\begin{bmatrix}U_{p}\\\ Y_{p}\\\
Y_{f_{L}}\end{bmatrix}g=\begin{bmatrix}u_{[t_{0},t_{0}+T_{p}-1]}\\\
y_{[t_{0},t_{0}+T_{p}-1]}\\\
y_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}+L-1]}\end{bmatrix}.$ (16)

###### Proof.

We first claim that (16) always has a solution. Since
$u_{[t_{0},t_{0}+T_{p}-1]}$ and $y_{[t_{0},t_{0}+T_{p}+T_{f}+L-1]}$ are
(different length) trajectory of (1), there exists a $T_{f}$ long sequence
$u^{*}$ that satisfies

$\begin{bmatrix}U_{p}\\\ Y_{p}\\\ U_{f}\\\
Y_{f_{L}}\end{bmatrix}g=\begin{bmatrix}u_{[t_{0},t_{0}+T_{p}-1]}\\\
y_{[t_{0},t_{0}+T_{p}-1]}\\\ u^{*}\\\
y_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}+L-1]}\end{bmatrix}$ (17)

by Lemmas 1 and 2, which proves the first claim.

We note that, for any solution $g$ to (16), the vector $U_{f}g$ is unique.
Indeed, if $g_{1}$ and $g_{2}$ are two solutions of (16), then
$g_{1}-g_{2}\in\mathcal{N}([U_{p}^{T},Y_{p}^{T},Y_{f_{L}}^{T}]^{T})$. This in
turn implies that $g_{1}-g_{2}\in\mathcal{N}(U_{f})$ by (12), so that
$U_{f}(g_{1}-g_{2})=0$.

Therefore, set $u_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}-1]}=U_{f}g$ for any solution
$g$ of (16). Then, it holds that

$\begin{bmatrix}U_{p}\\\ Y_{p}\\\ U_{f}\\\
Y_{f_{L}}\end{bmatrix}g=\begin{bmatrix}u_{[t_{0},t_{0}+T_{p}-1]}\\\
y_{[t_{0},t_{0}+T_{p}-1]}\\\ u_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}-1]}\\\
y_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}+L-1]}\end{bmatrix},$ (18)

which implies that the pair $u_{[t_{0},t_{0}+T_{p}+T_{f}-1]}$ and
$y_{[t_{0},t_{0}+T_{0}+T_{f}+L-1]}$ are trajectories of (1) by Lemma 2. ∎

Theorem 2 provides a foundation to answer the question posed in Introduction,
i.e., determining the input that generated a given output using previously
collected data, since it states that $T_{f}$ long input
$u_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}-1]}$ is uniquely determined for any pair
$u_{[t_{0},t_{0}+T_{p}-1]}$ and $y_{[t_{0},t_{0}+T_{p}+T_{f}+L-1]}$. Figure 1
illustrates this with color coded signals when $T_{p}=4$, $T_{f}=3$ and $L=2$.
Signals in black are required to determine the signal in red.

Figure 1: Graphical illustration of sequences in Theorem 2

The necessity of $T_{p}$ long pair of previous $u$ and $y$ is for the system
state $x(t_{0}+T_{p})$ identification. From this perspective, an
interpretation of the theorem is given as follows. Any vector $g$ belonging to

$\mathcal{G}\mathrel{\mathop{\mathchar
58\relax}}=\left\\{g\mathrel{\mathop{\mathchar
58\relax}}\begin{bmatrix}U_{p}\\\
Y_{p}\end{bmatrix}g=\begin{bmatrix}u_{[t_{0},t_{0}+T_{p}-1]}\\\
y_{[t_{0},t_{0}+T_{p}-1]}\end{bmatrix}\right\\}$

contains the information of the state $x(t_{0}+T_{p})$, so that the set

$\left\\{\begin{bmatrix}U_{f}\\\
Y_{f_{L}}\end{bmatrix}g\mathrel{\mathop{\mathchar
58\relax}}g\in\mathcal{G}\right\\}$

is a collection of all possible $T_{f}$-long input and $(T_{f}+L)$-long output
sequences starting from $x(t_{0}+T_{p})$. Thus, any observed
$y_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}+L-1]}$ must belong to this set, i.e., there
exists $g\in\mathcal{G}$ and

$y_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}+L-1]}=Y_{f_{L}}g.$

This $g$ is identified by solving (16). The corresponding
$u_{[t_{0}+T_{p},t_{0}+T_{p}+T_{f}-1]}$ is obviously given as in (15) and
unique.

The result of Theorem 2 is written in the form of algorithm. See Algorithm 1.

1: Setup:  Set $T_{p},T_{f},L\in\mathbb{Z}$ with $T_{p}\geq$ (observability
index of the system). Make Hankel matrices $U_{p}$, $U_{f}$, $Y_{p}$,
$Y_{f_{L}}$ from persistently exciting input of order $n+T_{p}+T_{f}+L$. Get
$u_{[-T_{p},-1]}$ and $y_{[-T_{p},-1]}$. Set $k=0$.

2: Input:  $y_{[k,k+T_{f}+L-1]}$

3: Solve $g$ for  $\displaystyle\begin{bmatrix}Y_{p}\\\ U_{p}\\\
Y_{f_{L}}\end{bmatrix}g=\begin{bmatrix}y_{[k-T_{p},k-1]}\\\
u_{[k-T_{p},k-1]}\\\ y_{[k,k+T_{f}+L-1]}\end{bmatrix}.$

4: Set $u_{[k,k+T_{f}-1]}=U_{f}g$.

5: Output:  $\hat{u}=u_{[k,k+T_{f}-1]}$

6: Update $k\leftarrow k+T_{f}$

7: Repeat from 2

Algorithm 1 Input Estimation

To initiate Algorithm 1, $u_{[-T_{p},-1]}$ and $y_{[-T_{p},-1]}$ must be
obtained. If the system is initially at rest, both can be set to zeros. An
example is provided.

######  Example 1.

As an example, we consider the system of (1) with

$\displaystyle A$ $\displaystyle=\begin{bmatrix}-0.3&0\\\
0&-0.5\end{bmatrix},\quad B=\begin{bmatrix}2&1\\\ -1&1\end{bmatrix},$ (19a)
$\displaystyle C$ $\displaystyle=\begin{bmatrix}1&2\\\ 1&0\end{bmatrix},\quad
D=\begin{bmatrix}0&0\\\ 0&0\end{bmatrix}.$ (19b)

The system has $L$-inverse with $L=1$. Applying $u_{[0,99]}$ as shown in
Figure 2(a) to the system with $x(0)=0$ yields $y_{[0,99]}$ as shown in Figure
2(b).

(a) Inputs

(b) Outputs

Figure 2: Inputs and outputs

To apply Algorithm 1, we set $T_{p}=2$, $T_{f}=3$, and $L=1$, and make Hankel
matrices with data generated by a 30 long random input sequence satisfying PE
of order $n+T_{p}+T_{f}+L=8$. Assuming the system initially at rest, Algorithm
1 begins iteration at time $k=3$ with $u_{[0,1]}=0_{4\times 1}$,
$y_{[0,1]}=0_{4\times 1}$ and the first input to the algorithm is $y_{[2,5]}$.
Denoting the outcome of Algorithm by $\hat{u}(k)$ at time $k$, we obtain
$\hat{u}_{[2,97]}$, which is shown in Figure 3 in red.

Figure 3: Estimation of inputs

Clear, $\hat{u}$ and $\hat{u}$ are identical with $L=1$ step delay.

Note that $\hat{u}$ is computed in batches of length $T_{f}=3$. Iteration $k$
gives $\hat{u}_{[3k,3k+2]}$.

Algorithm 1 receives $y_{[k,k+T_{f}+L-1]}$ as an input and computes
$\hat{u}=u_{[k,k+T_{f}-1]}$ at a time. This may be cumbersome to implement for
real-time applications. We provide Algorithm 2 tailored for real-time
application. Here, the input is modified to use $y$ available up to the
current instance and compute $\hat{u}$ of length 1.

1: Setup:  Set $T_{p},L\in\mathbb{Z}$ with $T_{p}\geq$ (observability index of
the system). Make Hankel matrices $U_{p}$, $U_{f}$, $Y_{p}$, $Y_{f_{L}}$ from
persistently exciting input of order $n+T_{p}+1+L$. Get $u_{[-T_{p}-L,-L-1]}$
and $y_{[-T_{p}-L,-1]}$. Set $k=0$.

2: Input:  $y(k)$

3: Solve $g$ for  $\displaystyle\begin{bmatrix}Y_{p}\\\ U_{p}\\\
Y_{f_{L}}\end{bmatrix}g=\begin{bmatrix}y_{[k-T_{p}-L,k-L-1]}\\\
u_{[k-T_{p}-L,k-L-1]}\\\ y_{[k-L,k]}\end{bmatrix}.$

4: Set $u(k-L)=U_{f}g$.

5: Output:  $\hat{u}(k)=u(k-L)$

6: Update $k\leftarrow k+1$

7: Repeat from 2

Algorithm 2 Input Estimation (with $T_{f}=1$)

######  Example 2.

Consider the system of (1) with

$\displaystyle A$ $\displaystyle=\begin{bmatrix}1&0.05&0&0.1&0.5&0\\\
0.05&1&0.05&0.05&0.1&0.05\\\ 0&0.05&1&0&0.05&0.1\\\
-0.2&0.1&0.05&0.8&0.1&0.05\\\ 0.1&-0.2&0.1&0.1&0.8&0.1\\\
0&0.1&-0.2&0.05&0.1&0.8\end{bmatrix},$ $\displaystyle B$
$\displaystyle=\begin{bmatrix}0.1&0\\\ 0.1&0\\\ 0&0.1\\\ 0.1&0.05\\\
-0.1&0.1\\\ 0&-0.1\end{bmatrix},\quad C^{T}=\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\\\ 0&0&0\\\ 0&0&0\\\ 0&0&0\end{bmatrix},$

and $D=0_{3\times 2}$. This system has $L_{0}$-delay inverse with $L_{0}=1$.

Algorithm 2 is used with $L=1$ and $T_{p}=2$. Figure 4 shows the input and
output sequence. Figure 5 shows the $\hat{u}$ along with $u$. Indeed they are
identical with one step delay.

(a) Inputs

(b) Outputs

Figure 4: Inputs and outputs Figure 5: Estimation of inputs

######  Example 3.

The proposed inversion of a linear system can also be used for finding a
constrained input that yields a desired output. Consider a desired output
$y^{*}_{[0,T_{f}+L-1]}$ and an input constraint set
$\mathcal{U}\subset\mathbb{R}^{m}$. Suppose that we want to find a constrained
input sequence $u_{[0,T_{f}-1]}$ at time $k=0$ when the previous records of
input $u_{[-T_{p},-1]}$ and output $y_{[-T_{p},-1]}$ are available. Then,
solve the minimization problem for $g$:

$\displaystyle\min_{g}\;$
$\displaystyle\|Y_{f_{L}}g-y^{*}_{[0,T_{f}+L-1]}\|^{2}$ subject to
$\displaystyle\begin{bmatrix}U_{p}\\\
Y_{p}\end{bmatrix}g=\begin{bmatrix}u_{[-T_{p},-1]}\\\
y_{[-T_{p},-1]}\end{bmatrix}$ $\displaystyle
U_{f}g\in\mathcal{U}\times\cdots\times\mathcal{U}$

and set $u_{[0,T_{f}-1]}=U_{f}g$.

##  IV Application to Disturbance Observer

Consider the block diagram of Figure 6, where $G(z)$ is the plant and the
block ‘DBINV’ is the data-based inversion block of $G(z)$ implemented by
Algorithm 2. The signals $u_{0}$, $y$, $d$, $u$ are, respectively, command
input, output, disturbance, and the actual input to the plant. Note that the
transfer function from $u_{0}$ to $y$ is

$Y(z)=G(z)[U_{0}(z)+(1-z^{-L})D(z)]$ (21)

where $Y$, $U_{0}$, and $D$ are $z$-transform of the corresponding signals,
which can also be seen from the figure by

$\displaystyle\hat{u}(k)=d(k-L)+\Delta(k-L),\quad\hat{d}(k)=d(k-L),$
$\displaystyle\Delta(k)=u_{0}(k)-\hat{d}(k)=u_{0}(k)-d(k-L).$

From (21) it is seen that, when $d(k)$ is slowly varying, the effect of $d$ is
approximately removed from the input. This is exactly the idea of the model-
based disturbance observer [10]. In this sense, we refer the proposed
structure as data-driven disturbance observer.

Figure 6: Data-driven Disturbance Observer

######  Example 4.

Consider the block diagram of Figure 6 where a realization of $G(z)$ is given
by the system of (19). We assume that the signal $u_{0}$ is the same as the
signal shown in Figure 2(a) and the signal $d$ is given in Figure 7(a). The
corresponding output is shown in Figure 7(b). Clearly, the effect of
disturbance is removed after some transients.

(a) Disturbances

(b) Outputs

Figure 7: Disturbances and outputs with disturbance observer

The simulation results agree with the derivation in (21) and we emphasize that
no information of $G(z)$ is directly used, and thus, the term of _data-driven
disturbance observer_ is justified.

##  V Conclusions

A data-based construction of inverse dynamics has been developed for discrete-
time LTI systems. The notion of $L$-delay inverse is invoked from the
literature, and it is combined with the system description method from
behavioral approach. The outcome is data-based representation of inverse
dynamics, which is similar to that of LTI systems, but differs in that input
is estimated with some delay. The result is applied to build disturbance
observers from no system information but collected input and output data that
satisfy a level of persistency of excitation condition. Applications and
extension this result seem to be immense.

## References

  * [1] J. C. Willems, P. Rapisarda, I. Markovsky, and B. L. M. De Moor, “A note on persistency of excitation,” Systems and Control Letters, vol. 54, no. 4, pp. 325–329, 2005.
  * [2] L. Huang, J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control for grid-connected power converters,” in the Proceedings of Conference on Decision and Control (CDC), pp. 8130-8135, 2019.
  * [3] J. Coulson, J. Lygeros, and F. Dörfler, “Data-enabled predictive control: In the shallows of the DeePC,” in the Proceedings of European Control Conference (ECC) pp. 307-312, 2019.
  * [4] P. G. Carlet, A. Favato, S. Bolognani, and F. Dörfler, “Data-driven predictive current control for synchronous motor drives,” in the Proceedings of Energy Conversion Congress and Exposition (ECCE), pp. 5148-5154, 2020.
  * [5] H. J. van Waarde, C. De Persis, M. K. Camlibel, and P. Tesi, “Willems’ fundamental lemma for state-space systems and its extension to multiple datasets,” IEEE Control Systems Letters, vol. 4, no. 3, pp. 602-607, 2020.
  * [6] C. De Persis and P. Tesi, “Formulas for data-driven control: Stabilization, optimality, and robustness,” IEEE Transactions on Automatic Control, vol. 65, no. 3, pp. 909–924, 2020.
  * [7] F. Dörfler, J. Coulson, and I. Markovsky, “Bridging direct & indirect data-driven control formulations via regularizations and relaxations,” IEEE Transactions on Automatic Control, 2022.
  * [8] I. Markovsky, J. C. Willems, S. Van Huffel, and B. De Moor, Exact and Approximate Modeling of Linear Systems: A Behavioral Approach. Philadelphia, PA, USA: SIAM, Mar. 2006.
  * [9] K. Ohishi, M. Nakao, K. Ohnishi, and K. Miyachi, “Microprocessor-controlled DC motor for load-insensitive position servo system,” IEEE Transactions on Industrial Electronics, vol. IE-34, no. 1, pp. 44–49, 1987.
  * [10] H. Shim, “Disturbance Observers,” In: Baillieul, J., Samad, T. (eds) Encyclopedia of Systems and Control, Springer, London, 2020. https://doi.org/10.1007/978-1-4471-5102-9_100068-1
  * [11] I. Markovsky and P. Rapisarda, “Data-driven simulation and control,” International Journal of Control, vol. 81, no. 12, pp. 1946-1959, 2008.
  * [12] M. Sain, J. Massey. “Invertibility of linear time-invariant dynamical systems.” IEEE Transactions on automatic control, vol. 14, no. 2, pp. 141-149, 1969.
  * [13] https://en.wikipedia.org/wiki/Rouche-Capelli_theorem