np.R
In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

### R code from vignette source 'np.Rnw'

###################################################
### code chunk number 1: np.Rnw:90-91
###################################################
options(prompt = "R> ", np.messages = FALSE, digits = 3, warn = -1)


###################################################
### code chunk number 2: np.Rnw:554-558
###################################################
library("np")
data("cps71")
model.par <- lm(logwage ~ age + I(age^2), data = cps71)
summary(model.par)


###################################################
### code chunk number 3: np.Rnw:570-576
###################################################
model.np <- npreg(logwage ~ age, 
                  regtype = "ll", 
                  bwmethod = "cv.aic",
                  gradients = TRUE,
                  data = cps71)
summary(model.np)


###################################################
### code chunk number 4: np.Rnw:592-593
###################################################
npsigtest(model.np)


###################################################
### code chunk number 5: np.Rnw:651-661 (eval = FALSE)
###################################################
## plot(cps71$age, cps71$logwage, xlab = "age", ylab = "log(wage)", cex=.1)
## lines(cps71$age, fitted(model.np), lty = 1, col = "blue")
## lines(cps71$age, fitted(model.par), lty = 2, col = " red")
## plot(model.np, plot.errors.method = "asymptotic")
## 
## plot(model.np, gradients = TRUE)
## lines(cps71$age, coef(model.par)[2]+2*cps71$age*coef(model.par)[3],
##       lty = 2,
##       col = "red")
## plot(model.np, gradients = TRUE, plot.errors.method = "asymptotic")


###################################################
### code chunk number 6: np.Rnw:666-691
###################################################
getOption("SweaveHooks")[["multifig"]]()
options(SweaveHooks = list(multifig = function() par(mfrow=c(2,2),mar=c(4,4,3,2))))

# Plot 1

plot(cps71$age,cps71$logwage,xlab="age",ylab="log(wage)",cex=.1)
lines(cps71$age,fitted(model.np),lty=1,col="blue")
lines(cps71$age,fitted(model.par),lty=2,col="red")

# Plot 2

plot(cps71$age,gradients(model.np),xlab="age",ylab="gradient",type="l",lty=1,col="blue")
lines(cps71$age,coef(model.par)[2]+2*cps71$age*coef(model.par)[3],lty=2,col="red")

# Plot 3

plot(cps71$age,fitted(model.np),xlab="age",ylab="log(wage)",ylim=c(min(fitted(model.np)-2*model.np$merr),max(fitted(model.np)+2*model.np$merr)),type="l")
lines(cps71$age,fitted(model.np)+2*model.np$merr,lty=2,col="red")
lines(cps71$age,fitted(model.np)-2*model.np$merr,lty=2,col="red")

# Plot 4

plot(cps71$age,gradients(model.np),xlab="age",ylab="gradient",ylim=c(min(gradients(model.np)-2*model.np$gerr),max(gradients(model.np)+2*model.np$gerr)),type="l",lty=1,col="blue")
lines(cps71$age,gradients(model.np)+2*model.np$gerr,lty=2,col="red")
lines(cps71$age,gradients(model.np)-2*model.np$gerr,lty=2,col="red")



###################################################
### code chunk number 7: np.Rnw:752-755
###################################################
cps.eval <- data.frame(age = seq(10,70, by=10))
predict(model.par, newdata = cps.eval)
predict(model.np, newdata = cps.eval)


###################################################
### code chunk number 8: np.Rnw:797-805
###################################################
data("wage1")
model.ols <- lm(lwage ~ female +
                married +
                educ +
                exper +
                tenure,
                data = wage1)
summary(model.ols)


###################################################
### code chunk number 9: np.Rnw:828-848
###################################################
model.ols <- lm(lwage ~ female +
                married +
                educ +
                exper +
                tenure,
                x = TRUE,
                y = TRUE,
                data = wage1)
X <- data.frame(wage1$female,
                wage1$married,
                wage1$educ,
                wage1$exper,
                wage1$tenure)
output <- npcmstest(model = model.ols, 
                    xdat = X, 
                    ydat = wage1$lwage, 
                    nmulti = 1,
                    tol = 0.1, 
                    ftol = 0.1)
summary(output)


###################################################
### code chunk number 10: np.Rnw:903-913
###################################################
#bw.all <- npregbw(formula = lwage ~ female +
#                  married +
#                  educ +
#                  exper +
#                  tenure,
#                  regtype = "ll",
#                  bwmethod = "cv.aic",
#                  data = wage1)
model.np <- npreg(bws = bw.all)
summary(model.np)


###################################################
### code chunk number 11: np.Rnw:946-981
###################################################
set.seed(123)
ii <- sample(seq(1, nrow(wage1)), replace=FALSE)
wage1.train <- wage1[ii[1:400],]
wage1.eval <- wage1[ii[401:nrow(wage1)],]
model.ols <- lm(lwage ~ female +
                married +
                educ +
                exper +
                tenure,
                data = wage1.train)
fit.ols <- predict(model.ols,
                   data = wage1.train,
                   newdata = wage1.eval)
pse.ols <- mean((wage1.eval$lwage - fit.ols)^2)
#bw.subset <- npregbw(formula = lwage ~ female +
#                     married +
#                     educ +
#                     exper +
#                     tenure,
#                     regtype = "ll",
#                     bwmethod = "cv.aic",
#                     data = wage1.train)
model.np <- npreg(bws = bw.subset)
fit.np <- predict(model.np,
                  data = wage1.train,
                  newdata = wage1.eval)
pse.np <- mean((wage1.eval$lwage - fit.np)^2)
bw.freq <- bw.subset
bw.freq$bw[1] <- 0
bw.freq$bw[2] <- 0
model.np.freq <- npreg(bws = bw.freq)
fit.np.freq <- predict(model.np.freq,
                       data = wage1.train,
                       newdata = wage1.eval)
pse.np.freq <- mean((wage1.eval$lwage - fit.np.freq)^2)


###################################################
### code chunk number 12: np.Rnw:1017-1020 (eval = FALSE)
###################################################
## plot(model.np,
##      plot.errors.method = "bootstrap",
##      plot.errors.boot.num = 25)


###################################################
### code chunk number 13: np.Rnw:1025-1028
###################################################
plot(model.np,
     plot.errors.method = "bootstrap",
     plot.errors.boot.num = 25)


###################################################
### code chunk number 14: np.Rnw:1073-1104
###################################################
data("birthwt", package = "MASS")
birthwt$low <- factor(birthwt$low)
birthwt$smoke <- factor(birthwt$smoke)
birthwt$race <- factor(birthwt$race)
birthwt$ht <- factor(birthwt$ht)
birthwt$ui <- factor(birthwt$ui)
birthwt$ftv <- factor(birthwt$ftv)
model.logit <- glm(low ~ smoke +
                   race +
                   ht +
                   ui +
                   ftv +
                   age +
                   lwt, 
                   family = binomial(link = logit),
                   data = birthwt)
model.np <- npconmode(low ~
                      smoke +
                      race +
                      ht +
                      ui +
                      ftv +
                      age +
                      lwt,
                      tol = 0.1,
                      ftol = 0.1,
                      data = birthwt)
cm <- table(birthwt$low, 
            ifelse(fitted(model.logit) > 0.5, 1, 0))
cm
summary(model.np)


###################################################
### code chunk number 15: np.Rnw:1132-1137
###################################################
data("faithful", package = "datasets")
f.faithful <- npudens(~ eruptions + waiting, data = faithful)
F.faithful <- npudist(~ eruptions + waiting, data = faithful)
summary(f.faithful)
summary(F.faithful)


###################################################
### code chunk number 16: np.Rnw:1142-1144 (eval = FALSE)
###################################################
## plot(f.faithful, xtrim = -0.2, view = "fixed", main = "")
## plot(F.faithful, xtrim = -0.2, view = "fixed", main = "")


###################################################
### code chunk number 17: np.Rnw:1149-1163
###################################################
getOption("SweaveHooks")[["multifig"]]()
options(SweaveHooks = list(multifig = function() par(mfrow=c(1,2),mar=rep(1.5,4))))
plot.output <- plot(f.faithful, xtrim=-0.2, view="fixed", main = "",plot.behavior="data")
# Retrieve data from plot() to create multiple figures
f <- matrix(plot.output$d1$dens, 50, 50)
plot.x1 <- unique(plot.output$d1$eval[,1])
plot.x2 <- unique(plot.output$d1$eval[,2])
persp(plot.x1,plot.x2,f,xlab="eruptions",ylab="waiting",zlab="Joint Density",col="lightblue", ticktype="detailed")

plot.output <- plot(F.faithful, xtrim = -0.2, view = "fixed", main="",plot.behavior="data")
# Retrieve data from plot() to create multiple figures
F <- matrix(plot.output$d1$dist, 50, 50)
plot.x1 <- unique(plot.output$d1$eval[,1])
plot.x2 <- unique(plot.output$d1$eval[,2])
persp(plot.x1,plot.x2,F,xlab="eruptions",ylab="waiting",zlab="Joint Distribution",col="lightblue", ticktype="detailed")


###################################################
### code chunk number 18: np.Rnw:1189-1200
###################################################
data("Italy")
fhat <- npcdens(gdp ~ year,
                tol = 0.1,
                ftol = 0.1,
                data = Italy)
summary(fhat)
Fhat <- npcdist(gdp ~ year,
                tol = 0.1,
                ftol = 0.1,
                data = Italy)
summary(Fhat)


###################################################
### code chunk number 19: np.Rnw:1207-1209 (eval = FALSE)
###################################################
## plot(fhat, view = "fixed", main = "", theta = 300, phi = 50)
## plot(Fhat, view = "fixed", main = "", theta = 300, phi = 50)


###################################################
### code chunk number 20: np.Rnw:1214-1229
###################################################
getOption("SweaveHooks")[["multifig"]]()
options(SweaveHooks = list(multifig = function() par(mfrow=c(1,2),mar=rep(1.25,4))))
plot.output <- plot(fhat, view="fixed", main="",plot.behavior="data")
# Retrieve data from plot() to create multiple figures
f <- matrix(plot.output$cd1$condens, 48, 50)
plot.y1 <- unique(plot.output$cd1$yeval)
plot.x1 <- unique(plot.output$cd1$xeval)
persp(as.integer(levels(plot.x1)),plot.y1,f,xlab="year",ylab="gdp",zlab="Conditional Density",col="lightblue", ticktype="detailed",theta=300,phi=50)

plot.output <- plot(Fhat, view="fixed", main="",plot.behavior="data")
# Retrieve data from plot() to create multiple figures
F <- matrix(plot.output$cd1$condist, 48, 50)
plot.y1 <- unique(plot.output$cd1$yeval)
plot.x1 <- unique(plot.output$cd1$xeval)
persp(as.numeric(plot.x1)+1951,plot.y1,F,xlab="year",ylab="gdp",zlab="Conditional Distribution",col="lightblue", ticktype="detailed",theta=300,phi=50)



###################################################
### code chunk number 21: np.Rnw:1259-1266
###################################################
bw <- npcdistbw(formula = gdp ~ year,
                tol = 0.1,
                ftol = 0.1,
                data = Italy)
model.q0.25 <- npqreg(bws = bw, tau = 0.25)
model.q0.50 <- npqreg(bws = bw, tau = 0.50)
model.q0.75 <- npqreg(bws = bw, tau = 0.75)


###################################################
### code chunk number 22: np.Rnw:1273-1280 (eval = FALSE)
###################################################
## plot(Italy$year, Italy$gdp, main = "", 
##      xlab = "Year", ylab = "GDP Quantiles")
## lines(Italy$year, model.q0.25$quantile, col = "red", lty = 1, lwd = 2)
## lines(Italy$year, model.q0.50$quantile, col = "blue", lty = 2, lwd = 2)
## lines(Italy$year, model.q0.75$quantile, col = "red", lty = 3, lwd = 2)
## legend(ordered(1951), 32, c("tau = 0.25", "tau = 0.50", "tau = 0.75"), 
##        lty = c(1, 2, 3), col = c("red", "blue", "red"))


###################################################
### code chunk number 23: np.Rnw:1289-1298
###################################################
plot(Italy$year, Italy$gdp, 
     main = "",
     xlab = "Year", 
     ylab = "GDP Quantiles")
lines(Italy$year, model.q0.25$quantile, col = "red", lty = 1, lwd = 2)
lines(Italy$year, model.q0.50$quantile, col = "blue", lty = 2, lwd = 2)
lines(Italy$year, model.q0.75$quantile, col = "red", lty = 3, lwd = 2)
legend(ordered(1951), 32, c("tau = 0.25", "tau = 0.50", "tau = 0.75"), 
       lty = c(1, 2, 3), col = c("red", "blue", "red"))


###################################################
### code chunk number 24: np.Rnw:1341-1347
###################################################
model.pl <- npplreg(lwage ~ female +
                    married +
                    educ +
                    tenure | exper,
                    data = wage1)
summary(model.pl)


###################################################
### code chunk number 25: np.Rnw:1367-1379
###################################################
model.index <- npindex(low ~
                       smoke +
                       race +
                       ht +
                       ui +
                       ftv +
                       age +
                       lwt,
                       method = "kleinspady",
                       gradients = TRUE,
                       data = birthwt)
summary(model.index)


###################################################
### code chunk number 26: np.Rnw:1399-1407
###################################################
model <- npindex(lwage ~ female +
                 married +
                 educ +
                 exper +
                 tenure,
                 data = wage1,
                 nmulti = 1)
summary(model)


###################################################
### code chunk number 27: np.Rnw:1433-1452
###################################################
model.ols <- lm(lwage ~ female +
                married +
                educ +
                exper +
                tenure,
                data = wage1)
wage1.augmented <- wage1
wage1.augmented$dfemale <- as.integer(wage1$female == "Male")
wage1.augmented$dmarried <- as.integer(wage1$married == "Notmarried")
model.scoef <- npscoef(lwage ~ dfemale +
                       dmarried +
                       educ +
                       exper +
                       tenure | female,
                       betas = TRUE,
                       data = wage1.augmented)
summary(model.scoef)
colMeans(coef(model.scoef))
coef(model.ols)


###################################################
### code chunk number 28: np.Rnw:1478-1480
###################################################
fit.lc <- npksum(txdat = cps71$age, tydat = cps71$logwage, bws = 2)$ksum/
  npksum(txdat = cps71$age, bws = 2)$ksum


###################################################
### code chunk number 29: np.Rnw:1485-1487 (eval = FALSE)
###################################################
## plot(cps71$age, cps71$logwage, xlab = "Age", ylab = "log(wage)")
## lines(cps71$age, fit.lc, col = "blue")


###################################################
### code chunk number 30: np.Rnw:1492-1494
###################################################
plot(cps71$age, cps71$logwage, xlab = "Age", ylab = "log(wage)", cex=.1)
lines(cps71$age, fit.lc, col = "blue")