Nothing
### 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")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.