if (FALSE) {
library(scam)
set.seed(2)
n <- 200
x1 <- runif(n)*6-3
f1 <- 3*exp(-x1^2) # unconstrained term
x2 <- runif(n)*4-1;
f2 <- exp(4*x2)/(1+exp(4*x2)) # monotone increasing smooth
f <- f1+f2
y <- f+rnorm(n)*0.2
dat <- data.frame(x1=x1,x2=x2,y=y)
b <- scam(y~s(x1,k=15,bs="cr")+s(x2,k=30,bs="mpi"),data=dat)
newd <- data.frame(x1=seq(-3,3,length.out=20),x2=seq(-1,3,length.out=20))
pred <- predict(b,newd)
pred
predict(b,newd,type="terms",se=TRUE)
## prediction on the original data...
pr <- predict(b,type="terms")
x<-sort(x1,index=TRUE)
old.par <- par(mfrow=c(2,2))
plot(x$x,(pr[,1])[x$ix],type="l",col=3,xlab="x1")
z<-sort(x2,index=TRUE)
plot(z$x,(pr[,2])[z$ix],type="l",col=3,xlab="x2")
plot(b,select=1,scale=0,se=FALSE)
plot(b,select=2,scale=0,se=FALSE)
par(old.par)
## linear extrapolation with predict.scam()...
set.seed(3)
n <- 100
x <- sort(runif(n)*3-1)
f <- exp(-1.3*x)
y <- rpois(n,exp(f))
dat <- data.frame(x=x,y=y)
b <- scam(y~s(x,k=15,bs="mpd"),family=poisson(link="log"),data=dat)
newd <- data.frame(x=c(2.3,2.7,3.2))
fe <- predict(b,newd,type="link",se=TRUE)
ylim<- c(min(y,exp(fe$fit)),max(y,exp(fe$fit)))
plot(c(x,newd[[1]]),c(y,NA,NA,NA),ylim=ylim,ylab="y",xlab="x")
lines(c(x,newd[[1]]),c(b$fitted,exp(fe$fit)),col=3)
## prediction on the original data...
pr <- predict(b)
plot(x,y)
lines(x,exp(pr),col=3)
## Gaussian model ....
## simulating data...
set.seed(2)
n <- 200
x <- sort(runif(n)*4-1)
f <- exp(4*x)/(1+exp(4*x)) # monotone increasing smooth
y <- f+rnorm(n)*0.1
dat <- data.frame(x=x,y=y)
b <- scam(y~ s(x,k=25,bs="mpi"),data=dat)
newd <- data.frame(x=c(3.2,3.3,3.6))
fe <- predict(b,newd)
plot(c(x,newd[[1]]),c(y,NA,NA,NA),ylab="y",xlab="x")
lines(c(x,newd[[1]]),c(b$fitted,fe),col=3)
### passing observed data + new data...
newd <- data.frame(x=c(x,3.2,3.3,3.6))
fe <- predict(b,newd,se=TRUE)
plot(newd[[1]],c(y,NA,NA,NA),ylab="y",xlab="x")
lines(newd[[1]],fe$fit,col=2)
lines(newd[[1]],fe$fit+2*fe$se.fit,col=3)
lines(newd[[1]],fe$fit-2*fe$se.fit,col=4)
## prediction with CI...
newd <- data.frame(x=seq(-1.2,3.5,length.out=100))
fe <- predict(b,newd,se=TRUE)
ylim<- c(min(y,fe$se.fit),max(y,fe$se.fit))
plot(newd[[1]],fe$fit,type="l",ylim=ylim,ylab="y",xlab="x")
lines(newd[[1]],fe$fit+2*fe$se.fit,lty=2)
lines(newd[[1]],fe$fit-2*fe$se.fit,lty=2)
## prediction on the original data...
pr <- predict(b)
plot(x,y)
lines(x,pr,col=3)
## bivariate example...
set.seed(2)
n <- 30
x1 <- sort(runif(n)); x2 <- sort(runif(n)*4-1)
f <- matrix(0,n,n)
for (i in 1:n) for (j in 1:n)
f[i,j] <- 2*sin(pi*x1[i]) +exp(4*x2[j])/(1+exp(4*x2[j]))
f <- as.vector(t(f));
y <- f+rnorm(length(f))*0.1
x11 <- matrix(0,n,n); x11[,1:n] <- x1; x11 <- as.vector(t(x11))
x22 <- rep(x2,n)
dat <- list(x1=x11,x2=x22,y=y)
b <- scam(y~s(x1,x2,k=c(10,10),bs="tesmi2"),data=dat,optimizer="efs")
par(mfrow=c(2,2),mar=c(4,4,2,2))
plot(b,se=TRUE); plot(b,pers=TRUE,theta = 80, phi = 40)
n.out <- 20
xp <- seq(0,1.4,length.out=n.out)
zp <- seq(-1,3.4,length.out=n.out)
xp1 <- matrix(0,n.out,n.out); xp1[,1:n.out] <- xp
xp1 <- as.vector(t(xp1)); xp2 <- rep(zp,n.out)
newd <- data.frame(x1=xp1,x2=xp2)
fe <- predict(b,newd)
fc <- t(matrix(fe,n.out,n.out))
persp(xp,zp,fc,expand= 0.85,ticktype = "simple",xlab="x1",
ylab="x2",zlab="f^",main="", theta = 80, phi = 40)
## obtaining a 'prediction matrix'...
newd <- data.frame(x1=c(-2,-1),x2=c(0,1))
Xp <- predict(b,newdata=newd,type="lpmatrix")
fv <- Xp%*% b$beta.t
fv
}
Run the code above in your browser using DataLab