McSpatial (version 1.1.1)

lwr: Locally Weighted Regression

Description

Estimates a model of the form y = f(x) using locally weighted regression. x can include either one or two variables. Returns estimated values, derivatives, and standard errors for both f(x) and df(x)/dx.

Usage

lwr(form,window=.25,bandwidth=0,kern="tcub",distance="Mahal",
  alldata=FALSE,data=NULL)

Arguments

form

Model formula

window

Window size. Default: 0.25.

bandwidth

Bandwidth. Default: not used.

kern

Kernel weighting function. Default is the tri-cube. Options include "rect", "tria", "epan", "bisq", "tcub", "trwt", and "gauss".

distance

Options: "Euclid", "Mahal", or "Latlong" for Euclidean, Mahalanobis, or "great-circle" geographic distance. May be abbreviated to the first letter but must be capitalized. Note: lwr looks for the first two letters to determine which vari

alldata

If alldata=T, each observation is used as a target value for x. When alldata=F, the function is estimated at a set of points chosen by the locfit program using an adaptive decision tree approach, and the akima<

data

A data frame containing the data. Default: use data in the current working directory.

Value

• target
  The target points for the original estimation of the function.
• ytarget
  The predicted values of y at the original target points.
• dtarget1
  The estimated derivatives dy/dx1 at the target points.
• dtarget2
  The estimated derivatives dy/dx2 at the target points. All zeros if the model has only one explanatory variable.
• ytarget.se
  Standard errors for the predicted values of y at the target points.
• dtarget1.se
  Standard errors for the derivatives dy/dx1 at the target points.
• dtarget2.se
  Standard errors for the derivatives dy/dx2 at the target points. All zeros if the model has only one explanatory variable.
• yhat
  The predicted values of y for the full data set.
• dhat1
  The estimated derivatives dy/dx1 for the full data set.
• dhat2
  The estimated derivatives dy/dx2 for the full data set. All zeros if the model has only one explanatory variable.
• yhat.se
  Standard errors for the predicted values of y for the full data set.
• dhat1.se
  Standard errors for the estimated derivatives dy/dx1 for the full data set.
• dhat2.se
  Standard errors for the estimated derivatives dy/dx2 for the full data set. All zeros if the model has only one explanatory variable.
• df1
  tr(L), a measure of the degrees of freedom used in estimation.
• df2
  tr(L'L), an alternative measure of the degrees of freedom used in estimation.
• sig2
  Estimated residual variance, sig2 = rss/(n-2*df1+df2).
• cv
  Cross-validation measure. cv = mean(((y-yhat)/(1-infl))^2), where yhat is vector of predicted values for y and infl is the vector of diagonal terms for L.
• gcv
  gcv = n*(n*sig2)/((n-nreg)^2), where sig2 is the estimated residual variance and nreg = 2*df1 - df2.
• infl
  A vector containing the diagonal elements of L.

Details

The estimated value of y at a target value $x_0$ is the predicted value from a weighted least squares regression of y on $x-x_0$ with weights given by $K(\psi/h)$, where $\psi$ is a measure of the distance between x and $x_0$ and h is the bandwidth or window. When x includes a single variable, $\psi = x-x_0.$ When x includes two variables, the method for specifying $\psi$ depends on the distance option. If distance="Mahal" or distance="Euclid", the ith row of the matrix X = (x1, x2) is transformed such that $x_i$ = sqrt(x_i * V * t(x_i)). Under the "Mahal" option, V is the inverse of cov(X). Under the "Euclid" option, V is the inverse of diag(cov(X)). By reducing x from two dimensions to one, this transformation leads again to the simple kernel weighting function $K((x- x_0 )/(sd(x)*h))$. The great circle formula is used to define K when distance = "Latlong"; in this case, the explanatory variable list must be specified as ~latitude+longitude (or ~lo+la or ~lat+long, etc), with the longitude and latitude variables expressed in degrees (e.g., -87.627800 and 41.881998 for one observation of longitude and latitude, respectively). The order in which latitude and longitude are listed does not matter and the function only looks for the first two letters to determine which variable is latitude and which is the longitude. It is important to note that the great circle distance measure is left in miles rather than being standardized. Thus, the window option should be specified when distance = "Latlong" or the bandwidth should be adjusted to account for the scale. The kernel weighting function becomes K(distance/h) under the "Latlong" option. h is specified by the bandwidth or window options. The intercept, $\alpha$, provides an estimate of y at $x_0$ and $\beta$ provides an estimate of the slope, dy/dx at $x_0$. When alldata=T, each data point in turn is used as a target point, $x_0$. Since each estimate is a linear function of all n values for y, the full set of estimates takes the form $\hat{y} = LY$, where L is an nxn matrix. Loader (1999) suggests two measures of the number of degrees of freedom used in estimation, df1 = tr(L) and df2 = tr(L'L), both of which are stored by lwr. The diagonal elements of tr(L) are stored in the array infl. Again following Loader (1999), the degrees of freedom correction used to estimate the error variance, $\hat{\sigma^2}$, is df = 2*df1 - df2. Let e represent the vector of residuals, $e = y-\hat{y}$. The estimated variance is $\hat{\sigma^2} = \sum_ie_i^2/(n-df)$. The covariance matrix is $$\hat{\sigma}^2(\sum_{i=1}^n Z_i K(\psi_i/h) Z_i^\top)^{-1}(\sum_{i=1}^n Z_i (K(\psi_i/h))^2 Z_i^\top )(\sum_{i=1}^n Z_i K(\psi_i/h) Z_i^\top)^{-1}.$$ where $Z = (1 \; x-x_0).$ Estimation can be very slow when alldata=T. When alldata=F, the package locfit is used to find a good set of target points at which to evaluate the function. See Loader (1999, section 12.2) for a description of the algorithm used to determine the target points. The akima package is then used to interpolate the coefficient estimates, the standard errors, and the values used to form df1 and df2. h can be specified to be either a fixed bandwidth or a window size set to a percentage of the sample size. Optionally, the lwr_grid command can be used to specify a vector of values for h with lwr picking the one that minimizes a criterion function. In general, the window option will be preferable because it provides more accurate estimates in regions where x is relatively sparse. Available kernel weighting functions include the following: lll{ Kernel Call abbreviation Kernel function K(z) Rectangular ``rect'' $\frac{1}{2} I(|z|

References

Cleveland, William S. and Susan J. Devlin, "Locally Weighted Regression: An Approach to Regression Analysis by Local Fitting," Journal of the American Statistical Association 83 (1988), 596-610. Loader, Clive. Local Regression and Likelihood. New York: Springer, 1999. McMillen, Daniel P., "Issues in Spatial Data Analysis," Journal of Regional Science 50 (2010), 119-141. McMillen, Daniel P., "Employment Densities, Spatial Autocorrelation, and Subcenters in Large Metropolitan Areas," Journal of Regional Science 44 (2004), 225-243. McMillen, Daniel P. and John F. McDonald, "A Nonparametric Analysis of Employment Density in a Polycentric City," Journal of Regional Science 37 (1997), 591-612. McMillen, Daniel P. and Christian Redfearn, "Estimation and Hypothesis Testing for Nonparametric Hedonic House Price Functions," Journal of Regional Science 50 (2010), 712-733. Pagan, Adrian and Aman Ullah. Nonparametric Econometrics. New York: Cambridge University Press, 1999. Silverman, A. W., Density Estimation for Statistics and Data Analysis, Chapman and Hall, New York (1986).

See Also

cparlwr cubespline fourier lwrgrid semip

Examples

# 1. Monte Carlo data
n = 1000
x <- runif(n,0,2*pi)
x <- sort(x)
ybase <- x - .1*(x^2) + sin(x) - cos(x) -.5*sin(2*x) + .5*cos(2*x)
sig = sd(ybase)/2
y <- ybase + rnorm(n,0,sig)
par(ask=TRUE)
plot(x,y)
lines(x,ybase,col="red")
fit <- lwr(y~x, window=.15)
# plot 95% confidence intervals for predicted y 
predse <- sqrt(fit$sig2 + fit$yhat.se^2)
lower <- fit$yhat + qnorm(.025)*predse
upper <- fit$yhat + qnorm(.975)*predse
plot(x, ybase, type="l", ylim=c(min(lower), max(upper)), 
  main="Estimated Function", xlab="x", ylab="y")
lines(x, fit$yhat, col="red")
lines(x, lower, lty="dashed", col="red")
lines(x, upper, lty="dashed", col="red")
legend("topleft", c("Base", "Predicted", "95 Percent CI"), 
 col=c("black", "red", "red"), lty=c("solid", "solid", "dashed"), lwd=1)

# plot 95%  confidence intervals for slopes
dxbase <- 1 - .2*x + cos(x) + sin(x) - cos(2*x) - sin(2*x)
lower <- fit$dhat1 + qnorm(.025)*fit$dhat1.se
upper <- fit$dhat1 + qnorm(.975)*fit$dhat1.se
plot(x, dxbase, type="l", ylim=c(min(lower), max(upper)), 
  main="Estimated Slopes", xlab="x", ylab="y")
lines(x, fit$dhat1, col="red")
lines(x, lower, lty="dashed", col="red")
lines(x, upper, lty="dashed", col="red")
legend("topright", c("Base", "Predicted", "95 Percent CI"), 
 col=c("black", "red", "red"),lty=c("solid", "solid", "dashed"), lwd=1)

# Derivative estimates with larger window size
fit <- lwr(y~x,window=.20)
lower <- fit$dhat1 + qnorm(.025)*fit$dhat1.se
upper <- fit$dhat1 + qnorm(.975)*fit$dhat1.se
plot(x, dxbase, type="l", ylim=c(min(lower), max(upper)), 
  main="Estimated Slopes", xlab="x", ylab="y")
lines(x, fit$dhat1, col="red")
lines(x, lower, lty="dashed", col="red")
lines(x, upper, lty="dashed", col="red")
legend("topright", c("Base", "Predicted", "95 Percent CI"), 
 col=c("black", "red", "red"), lty=c("solid", "solid", "dashed"), lwd=1)

#2. Population density data
library(maptools)
library(RColorBrewer)

cook <- readShapePoly(system.file("maps/CookCensusTracts.shp",
  package="McSpatial"))
cook$obs <- seq(1:nrow(cook))
# measure distance to Chicago city center
lmat <- coordinates(cook)
cook$LONGITUDE <- lmat[,1]
cook$LATITUDE  <- lmat[,2]
cook$DCBD <- geodistance(longvar=cook$LONGITUDE,latvar=cook$LATITUDE,
  lotarget=-87.627800,latarget=41.881998,dcoor=FALSE)$dist
# population density = population/acres,  acres = square mile x 640
cook$LNDENS <- log(cook$POPULATION/(cook$AREA*640))
densdata <- data.frame(cook[cook$POPULATION>0,])
par(ask=TRUE)

# lndens = f(longitude, latitude), weights are function of straight-line distance
fit <- lwr(LNDENS~LONGITUDE+LATITUDE,  window=.10, 
   distance="Latlong",data=densdata)
c(fit$df1, fit$df2, 2*fit$df1-fit$df2)
cook$lwrhat[densdata$obs] <- fit$yhat
brks <- seq(min(cook$lwrhat,na.rm=TRUE),max(cook$lwrhat,na.rm=TRUE),length=9)
spplot(cook,"lwrhat",at=brks,col.regions=rev(brewer.pal(9,"RdBu")),
   main="Log Density LWR Estimates")