\documentclass[11pt]{article}

\usepackage{fullpage}
\usepackage{psfig}
\usepackage{times}

\let\ul\underline
\let\bar\overline
\let\del\partial

\parindent=0pt
\parskip=0.5em

\begin{document}

\begin{center}
\Large\bf Nonlinear Regression: Some Ideas and a Little Guts
\end{center}

\ul{\em Nonlinear regression} models are models of the form
\[
        \begin{array}{rcl}
        y_i & = & f(\ul x_i ;~\ul\beta) + \epsilon_i; \\[0.5em]
        & & \epsilon_i  ~ \stackrel{iid}{\sim}  ~ N(0,\sigma^2), ~i=1,\ldots,n
        \end{array}
\]
where $\ul x_i = (x_{i1}, \ldots x_{iJ})$ as in linear regression,
and $\ul\beta=(beta_1, \beta_2, \ldots, \beta_{K})$ are parameters
that describe the function $f(\cdot;\ul\beta)$.  

Examples (all with $J=1$):
\begin{itemize}

\item Exponential model: $f(x;\alpha,\beta) = \displaystyle
\alpha e^{\beta x}$

\item Logistic model: $f(x;\alpha.\beta,\gamma) = \displaystyle
\frac{\alpha}{1 + \gamma e^{-\beta x}}$

\item Gompertz growth model: $f(x;\alpha,\beta,\gamma) = \displaystyle
\alpha e^{-\gamma e^{-\beta x}}$

\end{itemize}
etc. etc.

{\large\bf What is nonlinear regression doing in a linear regression course?}
\begin{itemize}
\item The nonlinear model 
\[
        y_i = f(x_{i1}, \ldots, x_{iJ)};~\ul\beta) + \epsilon_i, ~i=1,\ldots,n
\]
is a natural generalization of the linear model
\[
        y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_{K-1}
        x_{i(K-1)} + \epsilon_i, ~i=1,\ldots,n
\]
\item It is really linear regression in disguise!
\end{itemize}

{\large\bf Some Guts: How to estimate the $\ul\beta$'s.}

Want to find the values $\hat{\ul\beta} = (\hat\beta_1, \ldots,
\hat\beta_K)$ that minimize
\[
        RSS(\hat{\ul\beta}) = \sum_{i=1}^n (y_i - \hat y_i)^2 
        = \sum_{i=1}^n (y_i - f(\ul x_i ;~\hat{\ul\beta}))^2 
\]
This is a differential calculus problem again..  For each $\beta_j$, 
\[
        \frac{\del RSS(\ul\beta)}{\del\beta_j} 
        = -2\sum_{i=1}^n (y_i - f(\ul x_i ;~\ul\beta)) 
        \frac{\del f(\ul x_i; ~\ul\beta)}{\del\beta_j}
\]
Setting these equal to zero leads to the system of equations
\begin{eqnarray*}
        \sum_{i=1}^n \left(y_i - f(\ul x_i ;~\ul\beta)\right) 
        \frac{\del f(\ul x_i; ~\ul\beta)}{\del\beta_1} & = & 0 \\
        \sum_{i=1}^n \left(y_i - f(\ul x_i ;~\ul\beta)\right) 
        \frac{\del f(\ul x_i; ~\ul\beta)}{\del\beta_2} & = & 0 \\
        & \vdots &  \\
        \sum_{i=1}^n \left(y_i - f(\ul x_i ;~\ul\beta)\right) 
        \frac{\del f(\ul x_i; ~\ul\beta)}{\del\beta_K} & = & 0 
\end{eqnarray*}
We can arrange this in the form
$$
        F(X;\ul\beta)^T ~ \left[\ul y - \ul f(X;\ul\beta)\right] ~=~\ul 0
\eqno{(*)}
$$
where 
\[
        F(X;\ul\beta) = \left[
            \begin{array}{cccc}
              \frac{\del f(x_1;\ul\beta)}{\del\beta_1} & 
              \frac{\del f(x_1;\ul\beta)}{\del\beta_2} &
              \cdots & 
              \frac{\del f(x_1;\ul\beta)}{\del\beta_J} \\
              \frac{\del f(x_2;\ul\beta)}{\del\beta_1} & 
              \frac{\del f(x_2;\ul\beta)}{\del\beta_2} &
              \cdots & 
              \frac{\del f(x_2;\ul\beta)}{\del\beta_J} \\
              \vdots & \vdots & \ddots & \vdots \\
              \frac{\del f(x_n;\ul\beta)}{\del\beta_1} & 
              \frac{\del f(x_n;\ul\beta)}{\del\beta_2} &
              \cdots & 
              \frac{\del f(x_n;\ul\beta)}{\del\beta_J} \\
            \end{array}
            \right]
          ;~
        \ul y = \left[\begin{array}{c}
            y_1 \\ \vdots \\ y_n
            \end{array}\right]
          ;~ \mbox{and}~
        \ul f(X;\ul\beta) = \left[\begin{array}{c}
            f(\ul x_1;\ul\beta)\\ \vdots \\ f(\ul x_n;\ul\beta)
            \end{array}\right]
\]
Equation $(*)$ gives us the ``normal equations'', analogous to 
$$
        X^T(\ul y - X\ul\beta) ~=~ \ul 0
\eqno{(**)}
$$
in multiple regression.  However we can't solve $(*)$ in one operation
by pre-multiplying by $(X^TX)^{-1}$ to get
$\hat\beta=(X^TX)^{-1}X^T\ul y$ as in multiple regression.

Instead we solve ``iteratively'':  Suppose after $\ell$ steps we our
best guess for the parameters is $\ul\beta^{(\ell)}$, and we want to
improve it to $\ul\beta^{(\ell+1)}$.  We proceed as follows:
\begin{itemize}
\item Use Taylor's theorem to approximate $f(X;\ul\beta^{(\ell+1)})$ as
\[
        f(X;\ul\beta^{(\ell+1)}) \approx f(X;\ul\beta^{(\ell)}) 
        + F(X;\ul\beta^{(\ell)}) 
        \left(\ul\beta^{(\ell+1)} - \ul\beta^{(\ell)} \right)
\]
\item Rewrite $(*)$ as
\[
        F(X;\ul\beta^{(\ell)})^T ~ 
        \left[
          \left(\ul y - f(X;\ul\beta^{(\ell)}) \right)
        + F(X;\ul\beta^{(\ell)}) 
        \left(\ul\beta^{(\ell+1)} - \ul\beta^{(\ell)}\right)
        \right] ~\approx~\ul 0
\]
Note that this looks exactly like $(**)$, so that
$$
        \ul\beta^{(\ell+1)} - \ul\beta^{(\ell)} =
          \left(F(X;\ul\beta^{(\ell)})^TF(X;\ul\beta^{(\ell)})\right)^{-1}
          (F(X;\ul\beta^{(\ell)})^T\left(\ul y - f(X;\ul\beta^{(\ell)}) \right)
\eqno{(***)}
$$

\item Equivalently, we can obtain $\ul\beta^* = \ul\beta^{(\ell+1)} -
\ul\beta^{(\ell)}$ by solving the least-squares regression
problem
\[
        \ul y^* = X^*\ul\beta^* + \ul \epsilon,
\]
where $\ul y^* = \ul y - f(X;\ul\beta^{(\ell)})$ and $X^* =
F(X;\ul\beta^{(\ell)})$. 
\end{itemize}
The \ul{\em Gauss-Newton method} iterates $(***)$ until there is
essentially no difference between $\ul\beta^{(\ell)}$ and
$\ul\beta^{(\ell+1)}$.  

The Gauss-Newton method requires
\begin{itemize}
\item Derivatives of $f(\ul x_i;\ul\beta)$;
\item Good starting values; and 
\item Some tuning of the step sizes in $(***)$
\end{itemize}
to work in practice, but most methods of estimating nonlinear
regressions are based on variants of Gauss-Newton.

\clearpage

\vfill

{\large\bf Comparison of Some Linear and Nonlinear Regression Facts.}

\vfill

\hspace*{-0.5in}
{\def\arraystretch{3}\tiny\normalsize\let\ds\displaystyle
\begin{tabular}{|r|l|l|}
\hline
Feature   & Linear Regression & Nonlinear Regression 
\\
\hline
Model 
& $\ul y = X\ul\beta + \ul\epsilon$ 
& $\ul y = \ul f(X;\ul\beta) + \ul\epsilon$ 
\\
Normal equations
& $X^T(\ul y - X\ul\beta) ~=~ \ul 0$
& $F(X;\ul\beta)^T ~ \left[\ul y - \ul f(X;\ul\beta)\right] ~=~\ul 0$ 
\\
Solution 
& $\hat{\ul\beta} = (X^TX)^{-1}X^T\ul y$
& $\hat{\ul\beta} = \lim_{\ell\rightarrow\infty} \ul\beta^{(\ell)}$,
~~~~~where ~~~~ $\ul\beta^{(\ell+1)} ~ = $ \\[-1.5em]
& & $ \ul\beta^{(\ell)} ~ +~ 
\left(F(X;\ul\beta^{(\ell)})^TF(X;\ul\beta^{(\ell)})\right)^{-1}
(F(X;\ul\beta^{(\ell)})^T\left(\ul y - f(X;\ul\beta^{(\ell)}) \right)$
\\
Residual variance 
& $ s^2 = \ds\frac{\sum_i (y_i - \ul x_i\cdot \hat{\ul\beta})^2}{n-K}$
& $ s^2 = \ds\frac{\sum_i (y_i -  f(\ul x_i;\hat{\ul\beta}))^2}{n-K}$
\\
$\widehat{{\rm Var}(\hat{\ul\beta})}$ 
& $s^2 (X^TX)^{-1}$ 
& $s^2 \left(F(X;\hat{\ul\beta})^TF(X;\hat{\ul\beta})\right)^{-1}$
\\
$R^2$
& $1 - \ds\frac{\sum_i (y_i - \ul x_i\cdot \hat{\ul\beta})^2}
{\sum_i (y_i - \bar y)^2}$
& $1 - \ds\frac{\sum_i (y_i -  f(\ul x_i;\hat{\ul\beta}))^2}
{\sum_i (y_i - \bar y)^2}$
\\
$t$, $F$, SS's diagnostics
& Based directly on $X$, $\ul y$.
& Imitate linear regression with $X^*=F(X;\hat{\ul\beta})$, $\ul y$.
\\
\hline
\end{tabular}
}

\vfill

\clearpage

{\large\bf Example.}

The following table lists British census data on the percentage of
women with at least one child, arranged by women's age and birthyear.
The graph shows curves for each birthyear, that is, each curve
represents one column of the table (see RWG, pp. 168ff.).  Women with
the same birthyear are said to be in the same {\em birth cohort}.

\begin{center}
\begin{tabular}{crrrrrrr}
\hline
Women's & \multicolumn{7}{c}{Women's Birthyear} \\
\cline{2-8}
  Age &  1920 &  1930 &  1940 &  1945 &  1950 &  1955 &  1960 \\
\hline
  15  &    0  &    0  &    0  &    0  &    0  &    0  &    0 \\
  20  &    7  &    9  &   13  &   17  &   19  &   18  &   13 \\
  25  &   39  &   48  &   59  &   60  &   53  &   45  &   39 \\
  30  &   67  &   75  &   82  &   82  &   75  &   68  &   NA \\
  35  &   76  &   83  &   87  &   88  &   83  &   NA  &   NA \\
  40  &   78  &   86  &   89  &   90  &   NA  &   NA  &   NA \\
  45  &   NA  &   86  &   89  &   NA  &   NA  &   NA  &   NA \\
\hline
\hline
\end{tabular}

\vspace*{1.25in}
\mbox{\psfig{figure=childraw.eps,width=4in}}
\vspace*{-1.25in}

\end{center}

The SPLUS code that generated the plot above is given next.
The asymmetric growth over time in each curve suggests we fit the
curves using the Gompertz model.

\begin{verbatim}
> child _ read.table("child.dat",header=T)
> names(child)
[1] "Age"   "X1920" "X1930" "X1940" "X1945" "X1950" "X1955" "X1960"
> matplot(child$Age,as.matrix(child)[,-1],type="l")
> legend(locator(),legend=names(child)[-1],lty=1:7)
> title("Percentage of women with at least one child,\nby women's age \ 
and year of birth")
> title(xlab="Women's age",ylab="Percent")
\end{verbatim}

\clearpage

The Gompertz model looks like
\[
        y_i = f(x_i;\alpha,\beta,\gamma) + \epsilon_i
\]
where
\[
        f(x;\alpha,\beta,\gamma) 
        = \alpha e^{-\gamma e^{-\beta x_i}} .
\]
We can fit this to each birth cohort (that is, to the curve
representing women from the same birth year), using the SPLUS function
{\tt nls}, for ``nonlinear least squares''.  

Let us start with the 1940 birth cohort.  First we need decent starting
values.  Taking logs, rearranging, and taking logs again, we get
\begin{eqnarray*}
        \log f(x) & = & \log \alpha - \gamma e^{-\beta x} \\
        \beta x & = & \log \gamma - \log(\log(\alpha/f(x))
\end{eqnarray*}
Letting $x\rightarrow\infty$ in the first equation, we see that
$\alpha$ is the maximum value for $f(x)$, $\alpha\approx 90$ for the
1940 data.  If we take $x=0$ in the second equation to get rid of
$\beta$, we discover $\gamma$ should be quite large (since
$f(0)\approx0$ in the graphs); taking a cue from the
text\footnote{Since I don't have much experience with the Gompertz
  model, it would have been hard for me to see how large $\gamma$
  should have been, and it really seems to matter for this problem!},
we will take $\gamma=1000$.  Now put $x=30$, so $f(30)=82$ according
to the table, and look at the second equation again.  We get $\beta =
(\log(1000) - \log\log(90/82))/30 \approx 0.31$.

Now we are ready to estimate the first curve using SPLUS:
\begin{verbatim}
> nls1940 _ nls(X1940 ~ alpha * exp(-gamma*exp(-beta*Age)),
+ data=child,start=list(alpha=90,beta=0.31,gamma=1000),trace=T)
4.50236 : 90 0.31 1000 
0.4183 : 89.1071 0.309252 933.509 
0.412138 : 89.1039 0.30956 941.892 
0.412135 : 89.1039 0.309562 941.97 
> summary(nls1940)

Formula: X1940 ~ alpha * exp( - gamma * exp( - beta * Age))

Parameters:
           Value  Std. Error  t value 
alpha  89.103900  0.19578900 455.1010
 beta   0.309562  0.00358636  86.3165
gamma 941.970000 75.35320000  12.5007

Residual standard error: 0.320989 on 4 degrees of freedom

Correlation of Parameter Estimates:
       alpha   beta 
 beta -0.508       
gamma -0.464  0.993
> 1 - ( (0.320989)^2 * 4 ) / sum((child$X1940-mean(child$X1940))^2)
[1] 0.9999527  # = R^2 --- wow, that's pretty good!
>
> # Make some plots of the fit...
>
> plot(child$Age,child$X1940,xlab="",ylab="")
> x _ seq(15,45,length=50)
> lines(x,89.1*exp(-942*exp(-.31*x)))
> title("Percent of women born in 1940 with children, by age",
+ xlab="Age of women", ylab="Percent")
> plot(nls1940$fitted.values,child$X1940-nls1940$fitted.values,xlab="",ylab="")
> abline(h=0)
> title("Residuals from 1940 fit",xlab="Fitted values",ylab="Residuals")
> qqn(child$X1940-nls1940$fitted.values)
\end{verbatim}

\vspace*{1.25in}
\hfil\mbox{\psfig{figure=gompfit1940.eps,width=4in}}

\hfil\mbox{\psfig{figure=gompres1940.eps,width=3in}}
\mbox{\psfig{figure=gompqqn1940.eps,width=3in}}

\clearpage

\begin{center}
\Large\bf Breakout session
\end{center}

Here is the code for fitting the 1945 birth cohort data, using the
1940 parameter estimates as starting values.
\begin{verbatim}
> child45 _ na.omit(child[,c("Age","X1945")])
> nls1945 _ nls(X1945 ~ alpha * exp(-gamma*exp(-beta*Age)),
+ data=child45, start=list(alpha=89.1,beta=0.31,gamma=942), trace=T)
16.9186 : 89.1 0.31 942 
5.01028 : 89.7436 0.294916 590.214 
1.98967 : 90.398 0.280863 443.121 
0.121371 : 90.4283 0.281552 465.866 
0.118423 : 90.4253 0.281703 468.054 
0.118423 : 90.4253 0.281703 468.057 
> summary(nls1945)

Formula: X1945 ~ alpha * exp( - gamma * exp( - beta * Age))

Parameters:
           Value  Std. Error  t value 
alpha  90.425300  0.16066900 562.8050
 beta   0.281703  0.00221843 126.9830
gamma 468.057000 22.54640000  20.7597

Residual standard error: 0.198681 on 3 degrees of freedom

Correlation of Parameter Estimates:
       alpha   beta 
 beta -0.634       
gamma -0.587  0.993

> 1 - ( (??????)^2 * ?? ) / sum(child45$X1945-mean(child45$X1945)^2)
[1] ????????  # = R^2 
>
> # Make some plots of the fit...
>
> plot(child45$Age,child45$X1945,xlab="",ylab="")
> x _ seq(15,45,length=50)
> lines(x,??????)
> title("Percent of women born in 1945 with children, by age",
+ xlab="Age of women", ylab="Percent")
> plot(nls1945$fitted.values,child45$X1945-nls1945$fitted.values,
+ xlab="",ylab="")
> abline(h=0)
> title("Residuals from 1945 fit",xlab="Fitted values",ylab="Residuals")
> qqn(child45$X1945-nls1945$fitted.values)
\end{verbatim}

\clearpage

\begin{enumerate}

\item Explain {\em carefully} what {\tt child45} is, what it contains,
  and why it had to be made.

\vfill

\item Comment on the fitted nonlinear regression: How many iterations
  were required?  Are all the parameters estimated to be
  (significantly) different from zero?  Which parameter is most poorly
  estimated?  Are there any problems analogous to collinearity in this
  fit?

\vfill

\item If we had used the starting values {\tt start=list(alpha=1,
  beta=1, gamma=1)}, the {\tt nls()} function would have stopped with
  the message {\tt Problem in nls(X1945 \~{} alpha * exp( - gamma *
  exp( -..: singular gradient matrix}.  What is going wrong?

\vfill

\item The value for \verb|sum(child45$X1945-mean(child45$X1945)^2)|
is approximately 7529.  Fill in the ``?????'' and the ``??'' in the
  line that calculates $R^2$, and then calculate $R^2$.  Comment on
  this value.

\vfill

\clearpage

\item Fill in the ``??????'' in the {\tt lines()} command that makes
  the plot of the fitted Gompertz curve for the 145 birth cohort.


\vfill

\vspace*{1in}
\hfil\mbox{\psfig{figure=gompfit1945.eps,width=4in}}

\hfil\mbox{\psfig{figure=gompres1945.eps,width=3in}}
\mbox{\psfig{figure=gompqqn1945.eps,width=3in}}

\vspace*{-0.8in}

\item Do the plots look good?  What kind of plot could you make, to
  explore whether there are other local minima besides the ones {\tt
  nls()} found?

\vfill

\vfill

\vfill

\clearpage

\item The logistic function for nonlinear regression is 
\[
f(x;\alpha.\beta,\gamma) = \displaystyle
\frac{\alpha}{1 + e^{-\beta (x-\gamma)}}
\]
\begin{enumerate}
\item Do you expect the logistic fit to be better or worse than the
  Gompertz fit?  Why or why not?

\vfill

\item If you have time, write out the {\tt nls()} command you would
  use to fit the logistic function to the 1945 data (including
  starting values, if you can).  {\em Hint for starting values: }
\[
        \log\left(\frac{f(x)/\alpha}{1-f(x)/\alpha}\right) 
        =  {\beta (x - \gamma)}
\]


\vfill

\vfill

\vfill

\end{enumerate}
\end{enumerate}

\end{document}






