% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}


% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
\usepackage{float}
\usepackage[font=small, labelfont=bf,tableposition=top]{caption}
\DeclareCaptionLabelFormat{andfigure}{#1~#2 (L) and \figurename~\thefigure\ (R)}
%\usepackage{enumitem}
\usepackage{enumerate}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION  \else STUDENT VERSION \fi\\US Census Modeling}




\markboth{Population Modeling with Census Data}{Population Modeling with Census Data}

\author{Jean Marie Linhart \\
Central Washington University \\
Ellensburg WA 98926-7424 USA\\
jmlinhart@cwu.edu
}

\makeStitlePDFLaTex


\begin{abstract}
The United States Census, conducted every 10 years, gives data on the
United States population, that can be modeled with the exponential,
logistic, or Gompertz functions. Analysis of the per-unit population
growth leads to an understanding of the assumptions underlying each
model.  Parameters in the exponential and logistic models can be
estimated from per-unit population growth calculated from the data.
Parameters for all three models can also be estimated by nonlinear
least squares curve fitting, and these yield models with much higher
$R^2$ values, indicating better parameter estimates.  Models can be
evaluated by $R^2$ and on their predictive power.  Despite the
simplicity of the models, and the complicated nature of human
population changes, all of these models do a surprisingly good job of
matching the United States population data.
\end{abstract}


\section*{SCENARIO DESCRIPTION}
The United States conducts a census every 10 years, which gives us
data on the population from 1790 to the present
\cite{USCensusHist}, \cite{USPopData}, \cite{MestertonGibbons}. Table \ref{USAData} and Figure
\ref{USAData}, below, give the census data for the United States population.  Population
data for individual states and cities in the United States are also
available through the census, and data on the population of foreign
countries and cities are available on the world wide web from a
variety of sources. We are going to explore mathematical models for
populations, starting with the United States census data and following
up by exploring one other data set: a US state or city, a foreign
country or a foreign city, or possibly the population of the entire
world.  Choose your second population wisely; country borders are not
always static with time and so finding appropriate data can be
difficult.  If your interest is in Germany, whose borders were
rewritten several times before and at the conclusion World War II, you
will have an easier time with data for a city like Berlin, Stuttgart,
or Munich rather than trying to model the entire country. Each student
must choose a different population to model, in addition to the United
States population.

\begin{table}[H]
\begin{tabular}[b]{|c|r|}
\hline
Year & Population \\
\hline
1790 & 3929214 \\
1800 & 5308483 \\
1810 & 7239861 \\
1820 & 9638453 \\
1830 & 12866020 \\
1840 & 17069453 \\
1850 & 23191876 \\
1860 & 31443321 \\
1870 & 39818449 \\
1880 & 50155783 \\
1890 & 62947714 \\
1900 & 75994575 \\
1910 & 91972266 \\
1920 & 105710620 \\
1930 & 122775046 \\
1940 & 131669275 \\
1950 & 151325798 \\
1960 & 179323175 \\
1970 & 203302031 \\
1980 & 226545805 \\
1990 & 248718302 \\
2000 & 281424603 \\
2010 & 308745538 \\
\hline
\end{tabular}
\qquad\qquad
\includegraphics[width=4.5in]{1-066-USPopData.png}
\captionlistentry[figure]{United States Population Data}
\captionsetup{labelformat=andfigure}
\caption{United States Population Data \label{USAData}}
\end{table}

You will investigate applying several population growth models to your
data sets.  The models are


\begin{itemize}
\item Exponential model
$$\dfrac{dP}{dt} = rP \qquad \mbox{with solution } \qquad P(t) = P_0e^{rt}$$
\item Logistic model (for an accessible history see \cite{RefractModel})
$$\dfrac{dP}{dt} = r\left(1-\dfrac{P}{L}\right)P \qquad \mbox{with solution } \qquad P(t)  = \dfrac{L}{1 + \dfrac{(L- P_0)}{P_0}e^{-rt}}$$
\item Gompertz model
$$\dfrac{dP}{dt} = r\log\left(\dfrac{L}{P}\right)P \qquad\mbox{with solution }
\qquad P(t)  = {L}e^{-\log\left(\frac{L}{P_0}\right)e^{-rt}}$$
\end{itemize}

You will determine the best parameter values in each model to fit your
data.  The error used to evaluate the model fit to the data is the
$R^2$ value given by the {\bf coefficient of determination}, or the
{\bf adjusted} $R^2$ value.

Your modeling goals are:
\begin{enumerate}[I.]
\item Determine which model is theoretically the best fit for your  
population data based on the assumptions of the models.
\item Determine which model actually has the best fit overall for the data 
sets.  You cannot do this ``by eye''.  You must have a numerical
measure for best fit.
\item Using as much data as you can before 2010, find the best fit using 
only this data, and predict the population in 2010 (or another similar
year).  Determine which model(s) has (have) the best prediction(s).
The error associated with the prediction is the (signed) relative
error.
\item Determine if one model is clearly better than the others.
\item Determine what relevant factors are
omitted from the models.
\item Determine if there are issues with the data that affects the 
suitability of these models.
\end{enumerate}

Your goal is to write a project report which contains:
\begin{enumerate}[(a)]
\item A 150-250 word abstract summarizing your report.  The abstract should 
include the main conclusions of your report.
\item An introduction.  In the introduction you should:
\begin{enumerate}[i.]
\item Explain the theoretical differences between the models.  
To do this you will probably want to discuss growth rates, in
particular the per-unit population growth rate.  Make sure that you
define your terms carefully and use them consistently for all models.
\item Determine if any model or models make more theoretical sense than 
the others for the populations you are looking at. Predict which you
think will give the best fit and the reasons why.  If your reasoning
turns out to be incorrect, there is no penalty.  In fact, this is
good, as it reflects the scientific method at work.  We make hypotheses, and
sometimes they are not validated by the data.
\end{enumerate}
\item A methods section explaining how you will get parameters for the models, and how you will evaluate the models and their predictions.
\begin{enumerate}[i.]
\item If you use nonlinear least squares curve fitting, explain how you 
linearize the models and get initial values for the parameters.
Explain how you came up with your initial guess for $L$ in the
logistic and Gompertz models.
\item Explain how you evaluate the fits. 
\item Describe how you make and evaluate the predictions.
\end{enumerate}
\item A results and conclusions section for each data set in which you include:
\begin{enumerate}[i.]
\item One graph of all the data and the model fits.  A second graph 
of the data and the model fits made by omitting the most recent data
point, along with the model prediction and the most recent data
point.  You should be able to clearly read the title, axis labels, and
axis titles on your graphs in your report.
\item Tables with the parameter values you find associated with your model fits.  Example:

\medskip
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Fit & \multicolumn{3}{|c|}{US population (millions)}&\multicolumn{3}{|c|}{Reduced Model Predictions}\\
\hline
Model & $P_0$ & $L$ & $r$ & 2010 Data & Prediction & Error \%\\
\hline
Exponential & $\quad P_0 =  \quad$ &  & $r= $  & 308.7 & 341.4 & 10.6\\
\hline
Logistic & $P_0 = $ & $\quad L = \quad$ & $r =  $ & 308.7 &  & \\
\hline
Gompertz & $P_0= $& $L = $  & $r = $ & 308.7 & & \\
\hline
Fit & \multicolumn{3}{|c|}{Bryan TX population (thousands)}&\multicolumn{3}{|c|}{Reduced Model Predictions}\\
\hline
Exponential & $0.503$ &  & $0.0264$ & 65,660 & 76,394 & 16.4 \\
\hline
Logistic & $0.523$ & $10.17$ & $0.0421$ & 65,660 &  65,480 & -0.3\\
\hline
Gompertz & $0.512$ & $29.21$ & $0.0122$ & 65,660 & 68,036 & 3.6\\
\hline
\end{tabular}
\medskip

\item  Evaluate which models performed best for the data set.  Discuss 
any criticisms or comments about the population models based on your
data, knowledge of history, and/or ``uncommon sense''.  In other
words, you should write conclusions about your data.
\end{enumerate}
\item A general conclusions section for the entire project which incorporates:
\begin{enumerate}[i.]
\item Table(s) with
\begin{itemize}
\item The error associated with each model for each data set ($R^2$ values).
\item The signed relative error or percent signed relative error in predictions for each model and each data set.
\end{itemize}
Example:

\medskip
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
&\multicolumn{3}{|c|}{$R^2$} & \multicolumn{3}{|c|}{Prediction Error \%}\\
\hline
Data & Exponential  & Logistic & Gompertz & Exponential & Logistic & Gompertz\\
\hline
USA &  &  & & 10.6\%  &  &  \\
\hline
Bryan, TX & 0.950 & 0.998 & 0.998  & 16.35\% & -0.27\% & 3.62\% \\
\hline
\end{tabular}
\medskip

\item An evaluation of which model(s) performed best.  The table should allow for an easy comparison.
\item Analyze the meaning of the statement ``these models have too much memory to produce accurate predictions.'' 
\item Criticisms and potential improvements for these population models.  
What do these models fail to take into account?
\item Criticisms and potential improvements for the population data.  
Are there aspects of the data or data collection that impact the
appropriateness of the data for use with these models?
\end{enumerate}
\item A bibliography where you cite your sources of information.  Where did you get your data?  Where did you learn more about these models? 
\end{enumerate}


\section*{A STEP-BY-STEP GUIDE TO THE PROJECT}
What follows is information, exercises, and questions to help you think
about these models and prepare your report.


\subsection*{Picturing a good model}\label{PicModel}
 A graph of the data is provided in Figure \ref{DrawModel}.  We will
 explore some mathematical models for this data.  Questions and
 activities for you to complete are numbered below:

\begin{enumerate}[1)]
\item Draw a curve on the data that you think would be a reasonable model
for this data.
\item Draw a curve that you think would be a poor model
for this data.  
\item Using complete sentences, explain why you picked your reasonable
model and why you picked your poor model.
\newcounter{QuestionCounter}
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=4.5 in]{1-066-USPopData.png}
\vspace*{-1 cm}
\end{center}
\captionsetup{width=.8\linewidth}
\caption{Draw a curve on the data that you think would be a reasonable
model for the data.  Draw a curve on the data that you think would be
a poor model for the data.}\label{DrawModel}
\end{figure}

\subsection*{Possible Models}\label{PosModel}

Figure \ref{PopModels} contains graphs of several functions that might
be used to model populations. You are likely familiar with a linear
model and an exponential model.  The logistic and Gompertz models both
have two horizontal asymptotes, one at zero, and one at $L>0$.

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Judging by shape
characteristics alone, rank these models from best to worst for the United States data and explain your reasoning. Since the logistic and the Gompertz models have very similar shapes, you may put both in the same ranking.
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=4.5in]{1-066-PopModels1.png}

\vspace*{-1.5cm}
\end{center}
\captionsetup{width=.8\linewidth}
\caption{Possible models for the United States population data.\label{PopModels}}
\end{figure}

\subsection*{Introducing PPGR}\label{PPGRSec}

One way to evaluate potential competing models is to examine the
assumptions they give us about our data.  With a population model, we
often look at the Per (unit) Population Growth Rate (PPGR).  If the
population as a function of time is $P(t)$, then the PPGR is
\begin{equation}
\mbox{PPGR} = \dfrac{1}{P} \dfrac{dP}{dt}
\label{PPGR_eqn}
\end{equation} The PPGR measures how
much one unit of population, in our case, one human, contributes to
the growth (or reduction) in the population per unit time.


\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Analyze (\ref{PPGR_eqn}) and find the units on PPGR.
\item Say you are given
an estimate that human PPGR is equal to $\frac{1}{2}$ per year with
reasoning that women are half the population, and each
woman can give birth once a year, contributing a new member to the
population.  Do you think this PPGR is too high or too low? Why?
\item   Think
about the factors this estimate does not take into account, and
generate a better estimate (including units) for how big you
think the PPGR should be. What did you take into account? Explain how
you came up with your estimate.  
\item 
Should the PPGR be constant and unchanging, or are there situations or
scenarios in which you would expect it to increase or decrease?
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}



\subsection*{PPGR and modeling assumptions}\label{PPGRCalc}

Each of the possible models
given above can be represented as both a function and the differential
equation whose solution is that function.  Thus far we have considered
the following models

\begin{itemize}
\item Exponential model
$$\dfrac{dP}{dt} = rP \qquad \mbox{with solution } \qquad P(t) = P_0e^{rt}$$
\item Logistic model (for an accessible history see \cite{RefractModel})
$$\dfrac{dP}{dt} = r\left(1-\dfrac{P}{L}\right)P \qquad \mbox{with solution } \qquad P(t)  = \dfrac{L}{1 + \dfrac{(L- P_0)}{P_0}e^{-rt}}$$
\item Gompertz model
$$\dfrac{dP}{dt} = r\log\left(\dfrac{L}{P}\right)P \qquad\mbox{with solution }
\qquad P(t)  = {L}e^{-\log\left(\frac{L}{P_0}\right)e^{-rt}}$$
\end{itemize}

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Calculate the PPGR for each of these models. 
\item What happens to the PPGR in each of these models when there are very
few people?  If there aren't very many people, there should be plenty
of resources or raw materials in the environment to support human
life, so you can also think of this as a situation in which resources
are abundant or unlimited.  Explain what your result means for
humans in each model.  Do some of these limits make more sense for
human populations than others?
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

In both the Gompertz and in the logistic model, $L$ represents the
{\it limiting population}, meaning largest population the environment can
support, which also goes by the name {\it carrying capacity}.
In the graph of the logistic and
Gompertz equation, you can see that $P(t)$ grows to a horizontal
asymptote at $L$.

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Analyze the equations and determine the units on $L$.  
\item  What happens to the PPGR in the
logistic and Gompertz models as the population gets close to having
$L$ people?
\item  How does this relate to the behavior you see in the
graphs in Figure \ref{PopModels}? 
\item Does this make sense for human populations?
\item In the exponential and linear model, the population can keep getting
larger and larger. In both of these models what happens to the PPGR as
the population gets bigger and bigger?  Does this make sense for human
populations?
\item Based on your analysis of PPGR, rank the models from those that make
the most sense for human populations to those that make the least, and
explain the reasoning behind your rankings.
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

Both the shape and the modeling assumptions for the linear model make
it a bad match for our population data. Thus, we drop it from our discussion for the remainder of the project.


\subsection*{Solving ordinary differential equations}\label{Solving}

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Show how to use separation of variables to solve the differential
equation models to get the exponential, logistic and Gompertz modeling
functions.  Assume an initial condition $P(0) = P_0$.  You may also
need integration by substitution and/or partial fractions to integrate
the logistic and Gompertz differential equations.
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}
\subsection*{Parameter Estimation Using PPGR}\label{ParamEst1}

The United States population data \cite{USCensusHist},
\cite{USPopData}, \cite{USPopWiki} is given in Table \ref{USAData}. This
same data is also in the file {\tt USPopulation.csv}, which is a text
file tabular data in the form of``comma separated values'' that is
read by most spreadsheets and some computer algebra systems.  The
first parameter we need to do population modeling is $P_0 = P(0)$
which is the population at time zero.

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item What year should we pick to call 0?  We will refer to this as our zero
year.
\item  Why did you make this choice?
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate} 

Figure \ref{PPGR_P} shows a plot of a calculation for the PPGR
vs.~the population data.  Most of this discussion mimics the
discussion in \cite{MestertonGibbons}.  The furthest (blue) data point
to the right represents a PPGR of approximately 0.010 at a population
of $2.8\times10^8$.

\begin{figure}[h]
\begin{center}
\includegraphics[width=4.5in]{1-066-PPGR_P.png}
\vspace*{-0.5 cm}
\end{center}
\caption{A graph of PPGR vs $P$ for the United States population data.\label{PPGR_P}}
\end{figure}

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Describe a way to calculate PPGR from the data, and calculate a few
values by hand.  Can you think of different ways to calculate PPGR
from the data?
\item Use a computer algebra system to read in the United States
population data from {\tt USPopulation.csv} and use it to generate a
plot of PPGR vs.~$P$, like that in Figure \ref{PPGR_P}. Check
your mathematics by checking against the values you calculated by hand. Your
graph should look similar to that in Figure \ref{PPGR_P}, but may not
be exactly the same.
\item How can you use the
average PPGR and the best line approximation to the PPGR data to find
parameters for the exponential and logistic models? 
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

We can use least squares curve fitting to find the best fit line to
the PPGR vs.~$P$ data.  Let's start with the situation where you have
data of the form $(x_i, y_i)$, where $x$ is the independent variable,
and you want to model with a linear function $ y = f(x, a, b)$ that
depends on two parameters, $a$ (the slope of the line), and $b$ (the
$y$-intercept) as well as $x$.  Our problem is to minimize the sum of
the squared errors, also called the {\it sum of squared residuals}
(SS$_{\mbox{res}}$) $$\mbox{SS}_{\mbox{res}} =
\sum_i (y_i - f(x_i, a, b))^2$$ with respect to $a$ and $b$.  The
method of ordinary least squares uses linear algebra to minimize the
sum of the squared residuals.  Linear algebra methods for
least-squares minimization work not only for lines, but also for
functions $f$ that are polynomials in $x$; $f(x,a_0, \ldots, a_n) =
a_0 + a_1 x + a_2 x^2 +\cdots + a_n x^n$.  Ordinary least squares
returns unique values for the parameters ($a_0, \ldots, a_n$) that
minimize the sum of the squared residuals.


\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item  In Figure \ref{PPGR_P}, what is $x$, the independent variable? What are its units? What is $y$?  What are its units?
\item In Figure \ref{PPGR_P}, draw and label the distance $y_i - f(x_i, a, b)$ for a single data
point, where $f(x_i, a, b)$ is the function in represented by the  solid line (best line). 
\item What is the formula for the function $f(x, a, b)$ in terms of
$x$, $a$ and $b$?
\item Use your computer algebra system to find the best line to the 
$\mbox{PPGR}(P)$ vs.~$P$ data, as well as the average of PPGR$(P)$,
and plot both, as in Figure \ref{PPGR_P}.  There will be a
built-in routine to do this.
\item Find the values of the parameters for the exponential model, and 
those for the logistic model using your PPGR$(P)$ vs.~$P$ graph,
average, and best line.  For exponential and logistic, you need $P_0$
and $r$.  For logistic, you will also need $L$, the limiting
population.
\item Use your parameter values to make exponential and logistic
models for the data, and plot these with the USA population data.  Are
the models a good match for the data?
\item What would you expect a graph of PPGR($P$) vs $P$ to look like if 
the Gompertz model was best suited to the data?
\setcounter{enumi}{\theQuestionCounter}
\end{enumerate}

\subsection*{Nonlinear least squares parameter estimation}\label{ParamEst2}


There are many nonlinear functions, such as the exponential,
logistic, and Gompertz functions for which ordinary least squares, if
you can find a way to use it at all, does not give the parameter values
that minimize the sum of the squared residuals.  For these functions, we
turn to nonlinear least squares estimation.  Nonlinear least squares
estimation generally uses a stepping method to approximate the
minimum of the sum of the squared residuals and the associated parameter
values.  It may not succeed because there may not be a unique minimum
of the sum of the squared residuals, or the region around the minimum
might be very flat, making it difficult to find the location of the
minimum.

Similar to ordinary least squares, nonlinear least squares methods
also require you to provide your data ($x_i$'s and $y_i$'s), and the
function $ f(x, a, b, c)$, in which we assume there are three unknown
parameters to be estimated, $a$, $b$, and $c$.  Often, in order to get
a satisfactory answer, you must also provide a good set of starting
values for the parameters, so that the routine starts somewhere near
the minimum of the sum of squared residuals.


Fortunately, for our models, we can get a good set of starting values
for the parameters using ordinary least squares, but to do this we
have to rewrite the equations in the form of a line, $y = ax+b$,
where $a$ and $b$ are are the unknown parameters for slope and
$y$-intercept. This is easy for the exponential function, where taking
a logarithm of both sides transforms it to a linear equation
\begin{align*}
P(t)&= P_0 e^{rt},\\
\log(P(t))&= \log\left(P_0 e^{rt}\right),\\
\log(P(t))&= \log(P_0) +rt.
\end{align*}
The last equation is the linearization of the exponential equation. 

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item What are the parameters to be estimated in the exponential equation?
\item What do you use for your data if you want to use ordinary least squares estimation to find the parameters in the exponential model?  What are your $y_i$'s and $x_i$'s?
\item Using ordinary least squares estimation, you will get the slope,
$a$, and $y$-intercept, $b$.  What are $P_0$ and $r$ for the
exponential equation in terms of $a$ and $b$?
\item Use your computer algebra system to get values for $P_0$ and $r$ 
using ordinary least squares estimation.
\item Now use a nonlinear least-squares estimator to find improved values 
for $P_0$ and $r$.  Graph your data with the ordinary least squares
estimate for $P_0$ and $r$ and with the new values from the nonlinear
least-squares estimator.
\item With the United States data you should clearly see an
improvement in model fit from the original parameters found by
ordinary least squares and the parameters from nonlinear
least-squares.  Describe how the nonlinear least-squares model
appears to be better.
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

We can also linearize the logistic and the Gompertz equations.  Recall
that they are
\begin{equation}
P(t) = \dfrac{L}{1 + \frac{(L-P_0)}{P_0}e^{-rt}}\quad \mbox{(logistic)}\qquad \mbox{and} \qquad P(t) = Le^{-\log\left(\frac{L}{P_0}\right)e^{-rt}}\quad \mbox{(Gompertz).}
\end{equation}
For both, we first must look at the data and from the data come up
with a reasonable starting guess for the value of $L$.  Then, in
both equations, we solve for the part that is in the form $Ae^{-rt}$.
We know from our work with the exponential equation that we can take
logarithms to get a linearization.  In the logistic equation, the term
of the form $Ae^{-rt}$ is $\dfrac{(L-P_0)}{P_0}e^{-rt}$, and so our
logistic equation becomes in linearized form
\begin{align*}
\dfrac{L}{P(t)} -1 &= \dfrac{(L-P_0)}{P_0} e^{-rt},\\
\log\left(\dfrac{L}{P(t)} -1\right) &= \log\left(\dfrac{(L-P_0)}{P_0}\right) -rt.
\end{align*}

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Look at your data and choose a value for $L$ that is too small, one 
that is too big, and one that is reasonable. How do you know the value
that is too small is really too small?  How do you know the value that is too
big is really too big? Why is your reasonable choice for $L$ reasonable?
\item\label{nonlinearStart} Using your reasonable value for $L$, 
 what are the $y_i$ data values in ordinary least squares estimation
 with the linearized logistic equation?  What are the $x_i$ data values?
\item Once you have $a$ and $b$ from ordinary least squares with the 
linearization of the logistic equation, how do you get initial values
for $P_0$ and $r$ from these?
\item Use your computer algebra system to get values for $P_0$ and 
$r$ using ordinary least squares and your reasonable value for $L$.
\item\label{nonlinearEnd} Now use nonlinear least squares estimation 
to find improved values for $L$, $P_0$, and $r$.  Graph the data with
the logistic model created from the reasonable value for $L$, and the
ordinary least squares estimate for $P_0$ and $r$, and the logistic
model created from the nonlinear least squares estimation for $L$,
$P_0$, and $r$.  Which looks like the better model?  Why does it look
like the better model?
\item Now, find the linearization of the Gompertz model using 
your reasonable value for $L$, and repeat the work in questions
\ref{ParamEst2}.\ref{nonlinearStart} to \ref{ParamEst2}.\ref{nonlinearEnd}
for the Gompertz model.
\item Why do you think that nonlinear least squares gets us better
parameter values than the values for using the linearization and
ordinary least squares?  Hint: what does the logarithm do to large
data values and large errors in the linearized models?
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

\subsection*{Goodness of fit: coefficient of determination}\label{GoFit}

Sometimes we can see from a graph that one set of model parameters
gives us a better fit to data than another, but there are many cases
in which it is not clear which model gives a better fit to the data.
We prefer to use a quantitative measure for goodness of fit, such as
the {\it coefficient of determination} known as $R^2$.  The coefficient of
determination tells us how well a model does compared to a linear
model with an intercept only.  Our data is a set of $n$ values indexed
$(x_i, y_i)$.  Let $\bar{y} =
\frac{1}{n}\sum_i y_i$ be the mean of the dependent variable, and let
$f(x_i)$ be the values obtained by the model, then (\cite{WikiR2}, 
\cite{Stat2})
$$R^2 = 1-\dfrac{\sum_i (y_i - f(x_i))^2}{\sum_i (y_i - \bar{y})^2} =
1-\dfrac{SS_{\mbox{res}}}{SS_{\mbox{tot}}}$$

An $R^2$ of 1 indicates the model perfectly fits the data, and
generally, the closer $R^2$ is to 1, the better the model fit.  But
since $R^2$ is a comparison between a given model and a linear model
with an intercept only, you cannot directly compare two models with
their $R^2$ values \cite{R2Paper}.  Another problem is that increasing
the number of parameters in the model implicitly increases the power
of the model to explain the data, and so comparing the exponential
model which has two unknown parameters to the logistic or Gompertz
with three unknown parameters, is inherently unfair.  It is known that
$R^2$ will never decrease as parameters are included in a
regression model.

A related measure is the adjusted $R^2$ value which takes the number
of parameters, $p$, in the model compared to the number of data
points, $n$, into account.  Adjusted $R^2$ is denoted $\bar{R}^2$ and
a formula for it is \cite{WikiR2}, \cite{Stat2} $$\bar{R}^2 = 1 -
(1-R^2)\dfrac{n-1}{n-p-1} = R^2 - (1-R^2)\dfrac{p}{n-p-1}$$ where $n$
is the number of data points in the model, and $p$ is the number of
parameters in the model.
 
Interpretation of adjusted $R^2$ is similar to that for $R^2$, in that
values closer to 1 are considered better.  But adjusted $R^2$ also
cannot be used compare two models directly, since adjusted $R^2$ is a
measure of suitability comparing the given model with a linear model
with only an intercept \cite{R2Paper}.  One model might do slightly
better than a second in comparison to the linear model with only an
intercept, but that does not mean that it does better when directly
compared to the second model.

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Is $R^2$ or adjusted $R^2$ a better measure for comparing the
goodness of fit between the exponential, logistic and Gompertz models?
Why?
\item Calculate $R^2$ or  adjusted $R^2$ for the exponential, 
logistic and Gompertz models, and interpret your findings.  What does
this tell you about which model performs best?
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}


\subsection*{Predictions}\label{Predict}

Another way to evaluate how well the model performs is to investigate
how good it is at making predictions.  Since we do not want to wait
for the next census to be taken, this is evaluated by a ``Leave One
Out'' (LOO) method of omitting the final data point, creating a model
based on the remaining data (called the reduced data), then comparing
the prediction of the reduced/LOO model and the most recent data
point.  Comparison is usually done by signed relative error:
$$\mbox{Error} =
\dfrac{\mbox{predicted} -
\mbox{actual}}{\mbox{actual}}$$ 
The relative error always compares to the actual data value, not to
the model predicted value.

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item When is the signed relative error, as defined above, positive?  When is it negative?
\item Why do we calculate the signed relative error as
$$\mbox{ Error} = \dfrac{\mbox{predicted} - \mbox{actual}}{\mbox{actual}} \qquad
\mbox{instead of}\qquad \mbox{ Error} = \dfrac{\mbox{predicted} - \mbox{actual}}{\mbox{predicted}}$$
\item Why do we omit that final data point in finding the model parameters 
before evaluating the correctness of the prediction of the final data
point?  Why not just use the model based on all of the data?
\item Find model parameters using the reduced/LOO data set, leaving out the 
most recent data point.  Do you need to linearize again, or do you already have good starting values for the parameters?  What starting values for the parameters can you use?
\item Use the reduced/LOO model parameters to predict the most recent data 
point, and calculate the signed relative error.
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

A last way to evaluate a model's performance, at least for the
logistic and Gompertz models, is to look at the estimate or prediction
of $L$, the limiting population, and see if what you know about the
population gives you any insight as to which is more reasonable.

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item Do the predicted $L$ values from the logistic and Gompertz models 
give you any insight as to which is more reasonable?
\setcounter{QuestionCounter}{\theenumi}
\end{enumerate}

\subsection*{Conclusions}\label{conclusions}

\begin{enumerate}[1)]
\setcounter{enumi}{\theQuestionCounter}
\item With all the information and analysis available to you, which model
do you think performed the best, and why?  Is this the one you
predicted would be the best from your analysis of the model
hypotheses?
\item Someone said, ``These models have too much memory to produce accurate
predictions.''  Figure out what is meant by that and use it as part of
your criticism of these models.
\item What factors relevant to population size are omitted from these
population models?
\item Consider the population data sets.  Are there any aspects of the data
or the method of collecting the data that impact the appropriateness
of the data for modeling?
\end{enumerate}

\begin{thebibliography}{98}
\bibitem{R2Paper} Anderson-Sprecher, Richard. 1994.  Model Comparison and $R^2$. {\it The American Statistician}. 48:2, May 1994. pp. 113-117.

\bibitem{Stat2} Cannon, Ann R., et. al. 2013. {\it STAT2: Building Models for a World of Data}. New York: W. H. Freeman and Company.

\bibitem{LifeCourse} Cohan, C. L. and Cole, S. W. 2002. Life course transitions and natural disaster: marriage, birth and divorce following Hurricane Hugo. {\it Journal of Family Psychology}. 16:1, Mar 2002. pp. 14-25.

\bibitem{USCensusHist} Gibson, Campbell and Jung, Kay. 2002. {\it Population Division: Historical Census Statistics on Population Totals by Race 1790-1990...}. U.S. Census Bureau, Washington DC 20233. Available on the world wide web at \href{https://www.census.gov/content/dam/Census/library/working-papers/2002/demo/POP-twps0056.pdf}{https://www.census.gov/content/dam/Census/library/working-papers/2002/demo/POP-twps0056.pdf}. 

\bibitem{RefractModel} Kingsland, Sharon. 1982. The Refractory Model: The Logistic Curve and the History of Population Ecology.  {\it The Quarterly Review of Biology}. 57:1, Mar 1982. pp. 29-52.
\bibitem{MestertonGibbons} Mesterton-Gibbons, Mike. 2007. {\it A Concrete Approach to Mathematical Modeling}. Hoboken, NJ: John Wiley \& Sons, Inc.

\bibitem{USPopData} U. S. Census Bureau. 2010. {\it Resident Population Data (Text Version)}. U. S. Census Bureau, Washington DC 20233.  Available on the world wide web at \href{https://www.census.gov/2010census/data/apportionment-pop-text.php}{https://www.census.gov/2010census/data/apportionment-pop-text.php}.

\bibitem{WikiR2} Wikipedia Contributors. 2017. {\it Coefficient of Determination}. Wikipedia, The Free Encyclopedia, 31 August 2017, \href{https://en.wikipedia.org/wiki/Coefficient\_of\_determination}{https://en.wikipedia.org/wiki/Coefficient\_of\_determination}. Accessed 31 August 2017.

\bibitem{USPopWiki} Wikipedia Contributors. 2017. {\it Demography of the United States}. Wikipedia, The Free Encyclopedia, 24 August 2017, \href{https://en.wikipedia.org/wiki/Demography\_of\_the\_United\_States}{https://en.wikipedia.org/wiki/Demography\_of\_the\_United\_States}. Accessed 24 August 2017.


\end{thebibliography}


\teacher{
\section*{NOTES FOR TEACHERS}
This project was created for a senior-level mathematical modeling
course which had differential equations and linear algebra as
prerequisites.  Students may or may not be familiar with ordinary least
squares estimation from linear algebra.

The original intent of the project was to have students work in groups
of three to analyze three data sets, the United States population data
and two others, and to co-write a report.  One student would have
primary responsibility for writing the introduction.  A second student
would write the methods section and the results and conclusions for a
non-United States population data set.  The third student would write
results for the United States data and for the third non-United States
population data set, and a general conclusions section for the report.
Each student would be responsible for writing code to analyze the
United States data and create any plots or results required for the
report. The students not writing the introduction would also write
code to analyze another data set and create all relevant plots and
results required for the report.  While this is a group project that
requires group interaction, each student has individual responsibility
for certain sections, making it easy to determine individual grades.
The flow, consistency, and proofreading of the entire report, as well
as a group-authored abstract, is the basis for a shared portion of the
grade; yet it is usually clear from the individual portions where
adjustments up and down need to be made on the group grade.

As an alternative to writing a report, you might have students work in
groups or as a class to create presentations explaining various
portions of project or presenting data and results.

When doing group work I have students submit an individual group work
assessment in which they list their contributions to the project and
the contributions of the other group members.  They discuss what went
well or poorly, and how they might improve for the future.  Then they
are asked how they would distribute \$100 among the group as pay for
their work on the project and to explain why they would make such a
distribution.  This gives me further feedback on individual
contributions to the project, and gives students an opportunity to
reflect and learn from their experiences.

I also have students submit a short 250-500 word reflection on what they
learned from the project.

Even students well prepared for the class and for this project benefit
from structured thinking about the various aspects of it.  The
step-by-step guide is created to give the instructor concrete tasks to
assign and facilitate discussions with the students.  These may be
used piecewise by instructors wishing to explore only some aspects of
the project.

This project was created to give students experience using a nonlinear
least squares routine to perform parameter estimation, and with a
situation in which they are asked to select the ``best'' model when
the results they get will not be consistent.  There is no ``right
answer''; there are good arguments to be made to justify different
conclusions. There are many things that can go wrong along the way,
and each provides an opportunity for discussion and learning.

A related project on the logistic model and populations that also
covers many of these ideas in a different context is \cite{SIMIODE1}.

Python 2.7 code that goes through all of the steps in the step-by-step
guide and for all information and graphs required in the project
description is in the Jupyter Notebook file\\ {\tt
PopulationModelingWithCensusData.ipynb}.\footnote{You can download the Python 2.7 version of Anaconda (which contains SciPy and NumPy for scientific computing) at 
\url{https://www.anaconda.com/download/}.}  

\subsection*{Picturing a good model}
One mistake students make in thinking about modeling is to believe
that a model that exactly matches the data is a better model than one
that does not.  To this end, the instructor might want to provide
examples in which the model closely follows the data but has some
other undesired behavior.  An example, Figure \ref{PolyInterp}, shows
two least squares interpolating polynomials.  They are an excellent
fit to the data points, but have undesirable behavior outside the
years 1790-2010.

\begin{figure}[h]
\begin{center}
\includegraphics[width=4.5in]{1-066-USPop_polyfit.png}
\end{center}
\vspace*{-1cm}
\captionsetup{width=.8\linewidth}
\caption{Least squares (high-order) polynomial fits to the data. Both 
have undesirable behavior outside the years 1790-2010, although both
are excellent models for the data given.  Code to create this graph is in the file {\tt USADataPolynomials.ipynb}.
}\label{PolyInterp}
\end{figure}


\subsection*{Possible models}
A few features to contrast the logistic
and Gompertz model are visible here. The logistic model has a
symmetric approach to each of its two horizontal asymptotes.  The
Gompertz model, in contrast, has a more rapid increase from the lower
asymptote, and a slower approach to the upper asymptote.  This often
results in the Gompertz model predicting a higher value for $L$, the
upper asymptote, than the logistic model.

\subsection*{Introducing PPGR}
The PPGR suggested in item 3.1 is very high; some obvious criticisms
are that not all women are of reproductive age, not all women who are
able to will get pregnant, and some people will die.  Students can
compare their estimates with the actual PPGR from the United States population
data which they will compute in the next section.

PPGR varies, as, for instance, people tend to try to avoid reproducing
when their resources are very limited.  On the other hand, birthrates
are known to increase after a catastrophe \cite{LifeCourse}.

\subsection*{PPGR and modeling assumptions}
Each differential equation gives a formula for PPGR$(P)$. Students
should realize that they use the appropriate formula and take a limit
of PPGR as $P\to 0$ to find out what happens when there are very few
people.

The Gompertz and linear model have an infinite limit for PPGR as $P\to
0$, which could come about in a practical consequence either by the
gestational period getting faster and faster, or by humans giving
birth to more and more infants as $P$ gets closer to zero.  As both
options are preposterous, students should realize that this hypothesis
of the Gompertz and linear model is not appropriate for human populations.

The exponential model has a constant PPGR, which most students will
find unreasonable, preferring the behavior of PPGR$\to 0$ as $P\to L$.
Some, however, may make sense of this arguing that human ingenuity
will always find a way around any limiting resource.

\subsection*{Parameter estimation using PPGR}


Improvements are made from year-to-year in data collection, and, if
needed, estimation.  I would use this as an argument to use the most
recent population observation as $P_0$.  Using PPGR to do your
parameter estimation, you get a better exponential model and worse
logistic model using year 2010 (the most recent data point for the
United States data).

PPGR can be calculated using forward, backwards, or centered
differences.  If $\{t_i\}$ represents consecutive census years, with
corresponding populations represented as $\{P_i\}$, then a
forward-difference approximation for PPGR is $$\mbox{PPGR}(P_i) =
\dfrac{1}{P_i}\dfrac{P_{i+1} - P_i}{t_{i+1}-t_i}$$ This is the one
used in the graph of PPGR$(P)$ vs.~$P$ for the United States data in
Figure \ref{PPGR_P}.

The most frequent mistake students will make in getting their own
version of Figure \ref{PPGR_P} will be to instead plot PPGR$(P)$
vs.~year instead of vs.~$P$.

The exponential differential equation gives PPGR$(P) = r$, and so the
average of the PPGR values will give the estimate for $r$ in the
exponential model.  On the other hand, the logistic differential
equation gives PPGR$(P) = r\left(1-\dfrac{P}{L}\right)$, and so the
$y$-intercept of the best line will be $r$, and its $x$-intercept will
be $L$.

The models you get using the PPGR to calculate parameter values are not
particularly visually pleasing as a match for the data, see Figure
\ref{PPGRfit2010}. Using the year 1790 as the zero year for the United States data
sometimes improves model fits.  Using centered differences to
calculate PPGR always improves fits, but still does not give a good
match the United States data. Switching to nonlinear least squares parameter
estimation will make a substantial improvement in model fits.

\begin{figure}[h]
\begin{center}
\includegraphics[width=4.5 in]{1-066-USAPPGRfit_2010.png}
\end{center}
\vspace*{-1cm}
\captionsetup{width=.8\linewidth}
\caption{Model fits for United States data using PPGR calculated by forward 
differences. Year 2010 is used as the zero year.}
\label{PPGRfit2010}\end{figure}

If a Gompertz model was the most appropriate, we would expect
PPGR$(P)$ to increase as $P\to 0^+$.  While some values are clustered
above the best line in Figure \ref{PPGR_P}, you do not see an increase
in PPGR as $P$ decreases toward 0.

\subsection*{Nonlinear least squares parameter estimation}

A teacher should demonstrate or have students attempt using a
nonlinear estimator with the logistic or Gompertz model without
initial parameters, and make a graph of the result; it will be clear
that the estimation failed.  Python's scipy.optimize.curve\_fit
routine succeeded at finding parameters for the exponential model with
the United States data without initial guesses, but failed with the
logistic and Gompertz models.

The linearization of the Gompertz model, using a guessed
reasonable value for $L$ is
$$\log\left(\log\left(\dfrac{L}{P(t)}\right)\right) =
\log\left(\log\left(\dfrac{L}{P_0}\right)\right) - rt.$$ 
Students might try to use 
$$\log\left(\log\left(\dfrac{P(t)}{L}\right)\right) =
\log\left(-\log\left(\dfrac{L}{P_0}\right)\right) - rt,$$ 
which will result in interesting computer algebra results, since a
good value for $L$ is larger than the population values, making
$\frac{L}{P(t)}<1$, and so its logarithm will be negative. The second
logarithm will take the log of a negative number and most systems will
assign a NaN (not-a-number) for this value.

In each of the linearizations we use logarithms.  Logarithms have the
property of reducing large errors disproportionately, and so any
errors between the best fit approximation for the larger population
values effectively have less weight or significance than they do with
nonlinear least squares.

While using initial parameter values here is enough to get the
nonlinear parameter estimation working for the United States
population data, there are other data sets in which the nonlinear
estimator will fail to converge within the default maximum number of
iterations.  In particular, you should expect to have to readjust the
maximum number of iterations allowed for some data sets.  Overflow
warnings are also common.  These issues provide an opportunity for
discussion of how the estimators work; their strategy is a
systematic exploration of the parameter space, so that at each step
they move to new parameter values which come closer to the minimum.


\subsection*{Coefficient of determination}

Instead of giving seniors the detailed explanation of $R^2$ and
adjusted $R^2$ in the guide, consider having the class search the web
to read about $R^2$, the coefficient of determination, and have them
explain it to each other.  If necessary, pointed questions can lead
them to adjusted $R^2$ for use with models with differing numbers of
parameters.  Student volunteers can give preliminary explanations
of ideas on the board, with the instructor following up with
additional corrections and elaborations as required.

Values for $R^2$ are not directly comparable, yet this is often what
we do.  In this project, you will often find that you have two very
similar $R^2$ or adjusted $R^2$ values; for the United States data set the
adjusted $R^2$ values are $0.983$ for the exponential model, $0.997$
for logistic, and $0.999$ for Gompertz.  These are so close you really
cannot say that one model is definitively better than another.  All
three models are excellent choices for the data. We have to find
another way to distinguish which model is best.  However, if you find
that you have substantially different $R^2$ values, say $0.59$ and
$0.86$, you are justified in concluding that the model with the $R^2$
of $0.86$ is superior to the one with $R^2$ of $0.59$.

Rather than going through a long discussion on how to calculate and
interpret $R^2$ and adjusted $R^2$, you might point your students to a
routine to calculate the $R^2$ or adjusted $R^2$ and simply give them
instruction on how to interpret the results.

\subsection*{Predictions}

Here we see the Gompertz model edging out the logistic model for the
best prediction with the United States data set, although both are very close.
Gompertz has -1.3\% error, and logistic has -3.8\% error, as compared
to the exponential model with 10.6\% error.  Despite the unphysical
assumption of the Gompertz model, it has the best adjusted $R^2$ and the best
prediction (although not by a large margin) for the United States data.

The Gompertz model approaches its upper asymptote more slowly, and so
it will predict a higher $L$ value than logistic will.  This might be
a basis to argue that one model is better than another.  Logistic
predicts a limiting population of 490 million people for the United States,
whereas Gompertz predicts a limiting population of 1.4 billion, about
three times larger. The population of China is already 1.38 billion, with
land area slightly smaller than the United States.  This might be reason and
evidence to believe that Gompertz is more likely to be correct.

\subsection*{Conclusions}

The Gompertz model, based on $R^2$ and the predictions, edges out the
logistic model for the United States data, yet, analysis of PPGR makes it less
applicable to the data.  Generally logistic and Gompertz are in close
competition on population data sets for which performs better.
Students should be uncomfortable with this conclusion; they want a
clear winner and there is not one.

Students are asked to address the statement ``these models have too
much memory to produce accurate predictions,'' in the project
description.  The reason for this is that population changes largely
follow the Markovian property that ``the past influences the future
only through the present.''  In short, the population of 2010 is only
really influenced by the population data from 2000.  The population of
1900 has no effect.  Yet, in finding our parameters for these models,
we used the entire population history, and most of that history has no bearing on the future.  

At this stage, students should be thinking about what will change the
United States population from now to the next census.  The models take
proportional births and deaths into account, but do not consider
migration, wars, epidemics, or even birth control pills.
A truly realistic model would have to be very, very complicated.  The three
models we consider are very, very simple.


Likewise, there are issues with the United States data.  The 1790
census encompasses the European settlement to the continental US, and
does not consider the native American population. It wasn't until 1900
that the native American population was counted in the census. The
borders of the United States have changed substantially from 1790 to
the present, and so what exactly we mean by the population of the
United States has changed with them.

Another issue worth thinking about is the Three-Fifths Compromise, in which
African-American slaves were counted as 3/5 of a person for the
sake of apportionment, but careful examination of census reports,
e.g., \cite{Census1790} and \cite{Census1790Wiki}, indicates
this law was only used for apportionment and is not reflected in the
data totals.

It is amazing, given the limitations of these models and the
limitations in the data given to us, that these models perform as well
as they do.


\begin{thebibliography}{99}

\bibitem{Census1790} J. Phillips, et. al. 1793. {\it Return of the Whole Number of Persons within the Several Districts of the United States ...}. U.S. Census Bureau, Philadelphia, PA. Available on the world wide web at \href{https://www2.census.gov/prod2/decennial/documents/1790a.pdf}{https://www2.census.gov/prod2/decennial/documents/1790a.pdf}.  


\bibitem{Census1790Wiki} Wikipedia Contributors. 2017. {\it 1790 United States Census}. Wikipedia, The Free Encyclopedia, 11 September 2017, \href{https://en.wikipedia.org/wiki/1790\_United\_States\_Census}{https://en.wikipedia.org/wiki/1790\_United\_States\_Census}. 

\bibitem{SIMIODE1} Winkel, Brian. 2015. 1-018-T-LogisticPopGrowth.


\end{thebibliography}
}





\end{document}
