############
## scatterplot smoothing
x <- 1:1000/1000
mu <- exp(-400*(x-0.6)^2)+
5*exp(-500*(x-0.75)^2)/3+2*exp(-500*(x-0.9)^2)
y <- mu+0.5*rnorm(1000)
#fit with default knots
y.fit <- asp2(y~f(x,adap=TRUE))
plot(y.fit,residuals=TRUE,lwd=2,scb.lwd=2,scb.lty="dashed")
# with shaded confidence region.
# Use scb.lty="blank" to plot the shades only.
plot(y.fit,residuals=TRUE,shade=TRUE,scb.lty="blank")
if (FALSE) {
## Model with heteroscedastic errors
attach(mcycle)
y=accel
kn1 <- default.knots(times,20)
# fit model with constant residual variance
fit= asp2(accel~f(times,basis="os",degree=3,knots=kn1,adap=FALSE),
niter = 20, niter.var = 200)
# fit model with varying residual variance
fith=aspHetero(fit,times,tol=1e-8)
op <- par(mfrow = c(1,3))
plot(fit);plot(fith)
#sigma() returns the fitted varying residual variance
plot(sort(times),sigma(fith)[order(times)],type="l")
par(op)
## additive models
x1 <- 1:300/300
x2 <- runif(300)
mu1 <- exp(-400*(x1-0.6)^2)+
5*exp(-500*(x1-0.75)^2)/3+2*exp(-500*(x1-0.9)^2)
mu2 <- sin(2*pi*x2)
y2 <- mu1+mu2+0.3*rnorm(300)
y2.fit <- asp2(y2~f(x1,adap=TRUE)+f(x2,adap=TRUE))
# switch off adaptive fitting for the first function
y21.fit <- asp2(y2~f(x1,adap=FALSE)+f(x2,adap=TRUE))
op <- par(mfrow = c(2, 2))
plot(y2.fit)
plot(y21.fit)
par(op)
## scatterplot smoothing with specified knots and subknots
x <- 1:400/400
mu <- sqrt(x*(1-x))*sin((2*pi*(1+2^((9-4*6)/5)))/(x+2^((9-4*6)/5)))
y <- mu+0.2*rnorm(400)
kn <- default.knots(x,80)
kn.var <- default.knots(kn,20)
y.fit <- asp2(y~f(x,knots=kn))
y.fit2 <- asp2(y~f(x,knots=kn,var.knots=kn.var,adap=TRUE))
op <- par(mfrow = c(1, 2))
plot(y.fit)
plot(y.fit2)
par(op)
##################
#more examples
beta=function(l,m,x)
return(gamma(l+m)*(gamma(l)*gamma(m))^(-1)*x^(l-1)*(1-x)^(m-1))
f1 = function(x) return((0.6*beta(30,17,x)+0.4*beta(3,11,x))*1/0.958)
f2 = function(z) return((sin(2*pi*(z-0.5))^2)*1/.3535)
f3 = function(z)
return((exp(-400*(z-0.6)^2)+
5/3*exp(-500*(z-0.75)^2)+2*exp(-500*(z-0.9)^2))*1/0.549)
set.seed(1)
N <- 500
x1 = runif(N,0,1)
x2 = runif(N,0,1)
x3 = runif(N,0,1)
kn1 <- default.knots(x1,40)
kn2 <- default.knots(x2,40)
kn3 <- default.knots(x3,40)
kn.var3 <- default.knots(kn3,5)
y <- f1(x1)+f2(x2)+f3(x3)+0.3*rnorm(N)
# semiparametric model
fit1= asp2(y~x1+f(x2,basis="os",degree=3,knots=kn2,adap=FALSE)
+f(x3,basis="os",degree=3,
knots=kn3,var.knots=kn.var3,adap=FALSE),
niter = 20, niter.var = 200)
summary(fit1)
plot(fit1,pages=1)
# all effects flexible
# fit model with all smoothing parameters constant
fit2a= asp2(y~f(x1,basis="os",degree=3,knots=kn1,adap=FALSE)
+f(x2,basis="os",degree=3,knots=kn2,adap=FALSE)
+f(x3,basis="os",degree=3,knots=kn3,adap=FALSE),
niter = 20, niter.var = 200)
plot(fit2a,pages=1)
# fit model with last smoothing parameter adaptive
fit2b= asp2(y~f(x1,basis="os",degree=3,knots=kn1,adap=FALSE)
+f(x2,basis="os",degree=3,knots=kn2,adap=FALSE)
+f(x3,basis="os",degree=3,knots=kn3,adap=TRUE,
var.knots=kn.var3,var.basis="os",var.degree=3),
niter = 20, niter.var = 200)
# plot smoothing parameter function for covariate x3.
# Note that in the case of B-splines additional knots are added,
# see references.
plot(seq(0,1,length.out=42), fit2b$y.cov/fit2b$random.var[85:126],
ylab=expression(lambda(x3)),xlab="x3",type="l",lwd=3)
# compute 95% simultaneous confidence bands.
# You could skip this and use "fit2b" indstead of "scb2b" later on, however,
# if N is large, computing the SCBs various times can take some time
# if you don't need fitted values and bounds for all covariate points
# (can be computationally intensive due to large matrix dimensions),
# set calc.stdev=F such that these are not computed.
scb2b<- scbM(fit2b,calc.stdev=FALSE)
plot(scb2b,pages=1)
# plot only f(x2).
plot(scb2b,select=2,mfrow=c(1,1),lwd=3,ylab="f(x2)",xlab="x2")
# plot.scbm (and plot.asp) returns fitted values and confidence limits,
# if you only need the returned object set plot=FALSE
pscb=plot(scb2b,plot=FALSE)
# add pointwise confidence intervals to the plot
polygon(c(pscb$grid.x[[2]], rev(pscb$grid.x[[2]])),
c(pscb$fitted[[2]]+1.96*pscb$Stdev[[2]],
rev(pscb$fitted[[2]]-1.96*pscb$Stdev[[2]])),
col = grey(0.85), border = NA)
lines(pscb$grid.x[[2]],pscb$lcb[[2]],lty="dotted",lwd=3)
lines(pscb$grid.x[[2]],pscb$fitted[[2]],lwd=3)
lines(pscb$grid.x[[2]],pscb$ucb[[2]],lty="dotted",lwd=3)
# plot first derivative of f(x1).
# Useful to check statistical significance of certain features (such
# as bumps) in a curve.
scb2bdrv<- scbM(fit2b,drv=1,calc.stdev=FALSE)
plot(scb2bdrv,select=1)
#the following would give the same result
#x11();plot(fit2b,select=1,drv=1)
# different style
plot(scb2bdrv,select=1,scb.lty="blank",
shade=TRUE,shade.col="steelblue")
}
Run the code above in your browser using DataLab