\documentclass[12pt]{article}
% \VignetteIndexEntry{Cumulative Distribution Function (CDF) Deconvolution}
\author{Thomas Kincaid}
\title{Cumulative Distribution Function (CDF) Deconvolution}
\usepackage[colorlinks=true, citecolor=blue, urlcolor=red]{hyperref}
\usepackage{chicago}
\bibliographystyle{chicago}
\usepackage{Sweave}
\textwidth=6.5in
\textheight=9.0in
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.5in
\setlength{\parindent}{0in}
\setlength{\parskip}{.1in}

\begin{document}
\maketitle
\tableofcontents
\setkeys{Gin}{width=1.0\textwidth}

\section{Introduction}

This document presents deconvolution of a cumulative distribution function (CDF).  Convolution is a term that refers to a distribution that is a mixture of two or more component distributions.  Specifically, we are interested in the situation where the CDF is a mixture of the distribution for a variable of interest and measurement error.  The presence of measurement error will cause the CDF to occur across a wider range of variable values, which will lead to biased estimation of the CDF.  Deconvolution is the name of the procedure for removing the measurement error from the CDF.  In order to illustrate the deconvolution process, simulated data is used rather than data generated by a probability survey.  For additional discussion of measurement error and deconvolution see Kincaid and Olsen \citeyear{Kincaid+Olsen:2012}.

\section{Preliminaries}

The initial step is to use the library function to load the spsurvey package.  After the package is loaded, a message is printed to the R console indicating that the spsurvey package was loaded successfully.

Load the spsurvey package

<<preliminaries>>=
# Load the spsurvey package
library(spsurvey)

@
\begin{verbatim}
Version 3.5.0 of the spsurvey package was loaded successfully.

\end{verbatim}

\section{Load the Simulated Variables Data Set}
A data set was created that contains simulated species richness variables.  The initial variable was created using the Normal distribution with mean 20 and standard deviation 5.  Addition variables were created by adding 25\%, 50\%, and 100\% measurement error variance to the original variable, where the percent values refer to the variance of the original variable.  The data file also includes x-coordinates and y-coordinates, which were simulated using the Uniform distribution with a range from zero to one.  The data set contains 500 records.

The next step is to load the data set.  The data function is used to load the data set and assign it to a data frame named decon\_data.  The nrow function is used to determine the number of rows in the decon\_data data frame, and the resulting value is assigned to an object named nr. Finally, the initial six lines and the final six lines in the decon\_data data frame are printed using the head and tail functions, respectively.

Load the simulated variables data set

<<data>>=
# Load the data set and determine the number of rows in the data frame
data(decon_data)
nr <- nrow(decon_data)

@
Display the initial six lines in the data file.

<<data>>=
# Display the initial six lines in the data file
head(decon_data)

@
Display the final six lines in the data file.

<<data>>=
# Display the final six lines in the data file
tail(decon_data)

@
\section{Illustration of Extraneous Variance}

Measurement error is an example of extraneous variance, additional variance that is convoluted with the population frequency distribution for a variable of interest.  The impact of extraneous variance on the CDF will be investigated in this section.  As an initial step, the summary function is used to summarize the data structure of the species richness values for the original variable.

\begin{verbatim}
> summary(decon_data$Richness)
\end{verbatim}

<<summary, echo=FALSE>>=
# Use the summary function to summarize the data structure of the species
# richness variable
cat("\nSummarize the data structure of the species richness variable:\n")
summary(decon_data$Richness)

@
Note that the species richness values range between approximately 8 and 31.  By way of comparison, consider the range of value for the original variable plus 100\% additional measurement error.
  
\begin{verbatim}
> summary(decon_data$Richness_100)
\end{verbatim}

<<summary, echo=FALSE>>=
# Use the summary function to summarize the data structure of the species
# richness variable plus 100\% measurrement error
cat("\nSummarize the data structure of the species richness variable plus \n100% measurrement error:\n")
summary(decon_data$Richness_100)

@
The range of species richness values now has increased to 1 to 40, which reflects the bias that is introduced into the CDF.  To illustrate the bias induced by extraneous variance, CDFs for the species richness variables will be graphed.  The first step is to assign a vector of values at which the CDFs will be calculated.  The sequence (seq) function is used to create the set of values using the range 0 to 40.  Output from the seq function is assigned to an object named cdfvals.

\begin{verbatim}
> cdfvals <- seq(0,40,length=25)
\end{verbatim}

<<estimates, echo=FALSE>>=
# Assign a set of values at which to calculate CDF estimates
cdfvals <- seq(0,40,length=25)

@
The next step is to calculate a CDF estimate for the original species richness variable.  The cdf.est function in spsurvey will be used to calculate the estimate.  Since values for the species richness variables were not selected in a probability survey, the survey design weights argument (wgt) for the function is assigned a vector of ones using the repeat (rep) function.  Recall that the object named nr contains the number of rows in the decon\_data data frame.

\begin{verbatim}
> CDF_org <- cdf.est(z=decon_data$Richness,
+                    wgt=rep(1, nr),
+                    x=decon_data$xcoord,
+                    y=decon_data$ycoord,
+                    cdfval=cdfvals)
\end{verbatim}

<<estimates, echo=FALSE>>=
# Calculate a CDF estimate for the original species richness variable
CDF_org <- cdf.est(z=decon_data$Richness,
                   wgt=rep(1, nr),
                   x=decon_data$xcoord,
                   y=decon_data$ycoord,
                   cdfval=cdfvals)


@
Similarly, the cdf.est function will be used to calculate CDF estimates for the species richness variables that include measurement error.

\begin{verbatim}
> CDF_25 <- cdf.est(z=decon_data$Richness_25,
+                   wgt=rep(1, nrow(decon_data)),
+                   x=decon_data$xcoord,
+                   y=decon_data$ycoord,
+                   cdfval=cdfvals)

> CDF_50 <- cdf.est(z=decon_data$Richness_50,
+                   wgt=rep(1, nrow(decon_data)),
+                   x=decon_data$xcoord,
+                   y=decon_data$ycoord,
+                   cdfval=cdfvals)

> CDF_100 <- cdf.est(z=decon_data$Richness_100,
+                    wgt=rep(1, nrow(decon_data)),
+                    x=decon_data$xcoord,
+                    y=decon_data$ycoord,
+                    cdfval=cdfvals)
\end{verbatim}

<<estimates, echo=FALSE>>=
# Calculate a CDF estimate for the variable plus 25% extraneous variance
CDF_25 <- cdf.est(z=decon_data$Richness_25,
                  wgt=rep(1, nrow(decon_data)),
                  x=decon_data$xcoord,
                  y=decon_data$ycoord,
                  cdfval=cdfvals)

# Calculate a CDF estimate for the variable plus 50% extraneous variance
CDF_50 <- cdf.est(z=decon_data$Richness_50,
                  wgt=rep(1, nrow(decon_data)),
                  x=decon_data$xcoord,
                  y=decon_data$ycoord,
                  cdfval=cdfvals)

# Calculate a CDF estimate for the variable plus 100% extraneous variance
CDF_100 <- cdf.est(z=decon_data$Richness_100,
                   wgt=rep(1, nrow(decon_data)),
                   x=decon_data$xcoord,
                   y=decon_data$ycoord,
                   cdfval=cdfvals)

@
Density estimates facilitate visualizing the increased variance induced in a CDF by measurement error.  The ash1.wgt function in spsurvey, which implements the average shifted histogram algorithm \cite{Scott:1985}, will be used to calculate density estimates for each of the species richness variables.

\begin{verbatim}
> Density_org <- ash1.wgt(decon_data$Richness, nbin=25)
> Density_25 <- ash1.wgt(decon_data$Richness_25, nbin=25)
> Density_50 <- ash1.wgt(decon_data$Richness_50, nbin=25)
> Density_100 <- ash1.wgt(decon_data$Richness_100, nbin=25)

\end{verbatim}

<<estimates, echo=FALSE>>=
# Calculate density estimates
Density_org <- ash1.wgt(decon_data$Richness, nbin=25)
Density_25 <- ash1.wgt(decon_data$Richness_25, nbin=25)
Density_50 <- ash1.wgt(decon_data$Richness_50, nbin=25)
Density_100 <- ash1.wgt(decon_data$Richness_100, nbin=25)

@
CDF and density estimates are displayed in Figure~\ref{fig:extraneous}.  Estimates for the original species richness variable are displayed using red.  Estimates for variables containing measurement error are displayed using blue.  

\begin{figure}
\centering

<<figure1, fig=TRUE, echo=FALSE>>=
op <- par(mfrow=c(2,1), mgp=c(2.2,0.6,0), mar=c(3,3,2,0)+0.1)

plot(CDF_org$CDF$Value, CDF_org$CDF$Estimate.P, type="l", lwd=2, col="red",
     xlim=c(0, 40), ylim=c(0, 100), xaxt="n", xlab="", ylab="Percent")
title(main="CDF Estimates")
axis(side=1, at=c(0, 5, 10, 15, 20, 25, 30, 35, 40),
     labels=c("0", "5", "10", "15", "20", "25", "30", "35", "40"))
mtext("Species Richness", side=1, line=1.75, cex=par("cex"))
lines(CDF_25$CDF$Value, CDF_25$CDF$Estimate.P, lwd=2, col="blue")
lines(CDF_50$CDF$Value, CDF_50$CDF$Estimate.P, lty=5, lwd=2.5,
      col="blue")
lines(CDF_100$CDF$Value, CDF_100$CDF$Estimate.P, lty=4, lwd=2.5,
      col="blue")
legend(x="bottomright", inset=0.025, legend=c("Original Variable",
      "25% Added Variance", "50% Added Variance",
      "100% Added Variance"), lty=c(1,1,5,4), lwd=c(2,2,2.5, 2.5),
      bty="o", cex=1, col=c("red", "blue", "blue", "blue"))

plot(Density_org$x, Density_org$y, type="l", lwd=2, col="red",
     xlim=c(0, 40), ylim=c(0, 0.075), xaxt="n", xlab="",
     ylab="Resource Density")
title(main="Density Estimates")
axis(side=1, at=c(0, 5, 10, 15, 20, 25, 30, 35, 40),
     labels=c("0", "5", "10", "15", "20", "25", "30", "35", "40"))
mtext("Species Richness", side=1, line=1.75, cex=par("cex"))
lines(Density_25$x, Density_25$y, lwd=2, col="blue")
lines(Density_50$x, Density_50$y, lty=5, lwd=2.5, col="blue")
lines(Density_100$x, Density_100$y, lty=4, lwd=2.5, col="blue")

par(op)

@
\caption{CDF and Density Estimates for Species Richness Variables.}
\label{fig:extraneous}
\end{figure}
\newpage

\section{Deconvolution}

Deconvolution will be demonstrated using the species richness variable that includes 100\% additional measurement error variance.  The decon.est function in spsurvey will be used to calculate a deconvoluted CDF estimate, which uses an algorithm based on the procedure developed by Stefansky and Bay \citeyear{Stefansky+Bay:1996}.  Argument sigma for that function specifies the extraneous variance.  For an actual survey design, a value for sigma would be estimated from the survey data or obtained from some other source.  Since we are using simulated data, a direct estimate of sigma is available.  Specifically, the variance estimation function (var) will be used to estimate the additional variance for the species richness variable with added measurement error, and the result will be assigned to an object named extvar.

\begin{verbatim}
> extvar <- var(decon_data$Richness_100) - var(decon_data$Richness)
> CDF_decon <- cdf.decon(z=decon_data$Richness_100,
+                        wgt=rep(1,nr),
+                        sigma=extvar,
+                        x=decon_data$xcoord,
+                        y=decon_data$ycoord,
+                        cdfval=cdfvals)
\end{verbatim}

<<decon, echo=FALSE>>=
# Calculate the deconvoluted CDF estimate
extvar <- var(decon_data$Richness_100) - var(decon_data$Richness)
CDF_decon <- cdf.decon(z=decon_data$Richness_100,
                       wgt=rep(1,nr),
                       sigma=extvar,
                       x=decon_data$xcoord,
                       y=decon_data$ycoord,
                       cdfval=cdfvals)

@
The original CDF and deconvoluted CDF estimates for the species richness variable that includes 100\% additional measurement error variance are displayed in Figure~\ref{fig:deconvolution}.  The CDF estimate and confidence bounds for the original variable are displayed in red.  The deconvoluted CDF estimate and confidence bound are displayed in blue.  One consequence of deconvolution is increased confidence bound width, which can be observed in the CDFs in Figure~\ref{fig:deconvolution}.  For example, for a species richness value of 25, the confidence bounds for the original variable are 78.7\% to 84.5\%, which is a width of 5.8\%.  For the deconvoluted CDF, the confidence bounds are 74.2\% to 83.9\%, which is a width of 9.7\%.

\begin{figure}
\centering

<<figure2, fig=TRUE, echo=FALSE>>=
op <- par(mgp=c(2.2,0.6,0), mar=c(3,3,1,0)+0.1)

plot(CDF_100$CDF$Value, CDF_100$CDF$Estimate.P, type="l", lwd=2, col="red",
   xlim=c(0, 40), xaxt="n", xlab="Species Richness", ylab="Percent")
axis(side=1, at=c(0, 5, 10, 15, 20, 25, 30, 35, 40),
     labels=c("0", "5", "10", "15", "20", "25", "30", "35", "40"))
confcut <- 5
pctval <- c(confcut, 100-confcut)
tvalue <- CDF_100$CDF$Estimate.P >= pctval[1] &
          CDF_100$CDF$Estimate.P <= pctval[2]
x <-  interp.cdf(pctval, CDF_100$CDF$Estimate.P, CDF_100$CDF$Value)
ylow <- interp.cdf(pctval, CDF_100$CDF$Estimate.P, CDF_100$CDF$LCB95Pct.P)
yhi <- interp.cdf(pctval, CDF_100$CDF$Estimate.P, CDF_100$CDF$UCB95Pct.P)
value <- c(x[1], CDF_100$CDF$Value[tvalue], x[2])
lower <- c(ylow[1], CDF_100$CDF$LCB95Pct.P[tvalue], ylow[2])
upper <- c(yhi[1], CDF_100$CDF$UCB95Pct.P[tvalue], yhi[2])
lines(value, lower, lty=3, lwd=2.5, col="red")
lines(value, upper, lty=3, lwd=2.5, col="red")
lines(CDF_decon$CDF$Value, CDF_decon$CDF$Estimate.P, lwd=2, col="blue")
tvalue <- CDF_decon$CDF$Estimate.P >= pctval[1] &
          CDF_decon$CDF$Estimate.P <= pctval[2]
x <-  interp.cdf(pctval, CDF_decon$CDF$Estimate.P, CDF_decon$CDF$Value)
ylow <- interp.cdf(pctval, CDF_decon$CDF$Estimate.P, CDF_decon$CDF$LCB95Pct.P)
yhi <- interp.cdf(pctval, CDF_decon$CDF$Estimate.P, CDF_decon$CDF$UCB95Pct.P)
value <- c(x[1], CDF_decon$CDF$Value[tvalue], x[2])
lower <- c(ylow[1], CDF_decon$CDF$LCB95Pct.P[tvalue], ylow[2])
upper <- c(yhi[1], CDF_decon$CDF$UCB95Pct.P[tvalue], yhi[2])
lines(value, lower, lty=3, lwd=2.5, col="blue")
lines(value, upper, lty=3, lwd=2.5, col="blue")
legend(x="bottomright", inset=0.05, legend=c("Original CDF",
       "Confidence Bounds", "Deconvoluted CDF", "Confidence Bounds"),
       lty=c(1,3,1,3), lwd=c(2,2.5,2,2.5), bty="o", cex=1,
       col=c("red", "red", "blue", "blue"))

par(op)

@
\caption{Original and Deconvoluted CDF Estimates for a Species Richness Variable with Added Measurement Error.}
\label{fig:deconvolution}
\end{figure}

\newpage
\bibliography{spsurvey}

\end{document}
