##=======
## Gaussian model, two smooth terms: unconstrained and increasing...
## simulating data...
require(scam)
set.seed(4)
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
y <- f1+f2 +rnorm(n)*.5
dat <- data.frame(x1=x1,x2=x2,y=y)
## fit model, get results, and plot...
b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi"),data=dat)
b
summary(b)
plot(b,pages=1,shade=TRUE)
##=======
## Gaussian model, two smooth terms: increasing and mixed (decreasing and convex)...
## simulating data...
set.seed(4)
n <- 200
x1 <- runif(n)*4-1;
f1 <- exp(4*x1)/(1+exp(4*x1)) # increasing smooth
x2 <- runif(n)*3-1;
f2 <- exp(-3*x2)/15 # decreasing and convex smooth
y <- f1+f2 + rnorm(n)*.4
dat <- data.frame(x1=x1,x2=x2,y=y)
## fit model, results, and plot...
b <- scam(y~ s(x1,bs="mpi")+s(x2, bs="mdcx"),data=dat)
summary(b)
plot(b,pages=1,scale=0,shade=TRUE)
##=======
if (FALSE) {
## using the extended Fellner-Schall method for smoothing parameter selection...
b0 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="efs")
summary(b0)
## using the extended Fellner-Schall method for smoothing parameter selection,
## and BFGS for model coefficient estimation...
b0 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer=c("efs","bfgs"))
summary(b0)
## using optim() for smoothing parameter selection...
b1 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="optim")
summary(b1)
b2 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="optim",
optim.method=c("BFGS","fd"))
summary(b2)
## using nlm()...
b3 <- scam(y~ s(x1,bs="mpi")+s(x2,bs="mdcx"),data=dat,optimizer="nlm")
summary(b3)
}
##=======
## Poisson model ....
## simulating data...
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 <- rpois(n,exp(f))
dat <- data.frame(x1=x1,x2=x2,y=y)
## fit model, get results, and plot...
b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi"),
family=poisson(link="log"),data=dat,optimizer=c("efs","bfgs"))
summary(b)
plot(b,pages=1,shade=TRUE)
scam.check(b)
## Gamma model...
## simulating data...
set.seed(6)
n <- 300
x1 <- runif(n)*6-3
f1 <- 1.5*sin(1.5*x1) # unconstrained term
x2 <- runif(n)*4-1;
f2 <- 1.5/(1+exp(-10*(x2+.75)))+1.5/(1+exp(-5*(x2-.75))) # increasing smooth
x3 <- runif(n)*6-3;
f3 <- 3*exp(-x3^2) # unconstrained term
f <- f1+f2+f3
y <- rgamma(n,shape=1,scale=exp(f))
dat <- data.frame(x1=x1,x2=x2,x3=x3,y=y)
## fit model, get results, and plot...
b <- scam(y~s(x1,bs="ps")+s(x2,k=15,bs="mpi")+s(x3,bs="ps"),
family=Gamma(link="log"),data=dat,optimizer=c("efs","bfgs"))
b
summary(b)
old.par <- par(mfrow=c(2,2))
plot(b,shade=TRUE)
par(old.par)
## run with unpenalized terms...
b1 <- scam(y~s(x1,bs="ps")+s(x2,k=15,bs="mpi",fx=TRUE)+s(x3,bs="ps",fx=TRUE),
family=Gamma(link="log"),data=dat,optimizer=c("efs","bfgs"))
summary(b1)
if (FALSE) {
## bivariate example...
simu <- function(x,z) -exp(4*x)/(1+exp(4*x))+2*sin(pi*z)
xs <- seq(-1,3,length=30); zs <- seq(0,1,length=30)
pr <- data.frame(x=rep(xs,30),z=rep(zs,rep(30,30)))
truth <- matrix(simu(pr$x,pr$z),30,30)
set.seed(2)
n <- 500
x <- runif(n)*4-1
z <- runif(n)
f <- simu(x,z)
y <- f + rnorm(n)*.2
## fit model ...
b <- scam(y~s(x,z,k=c(10,10),bs=c("tesmd1","ps")),optimizer="efs")
summary(b)
old.par <- par(mfrow=c(2,2),mar=c(4,4,2,2))
plot(b,se=TRUE)
plot(b,pers=TRUE,theta = 30, phi = 40);title("tesmd1")
plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data",pch=".",cex=3)
persp(xs,zs,truth,theta = 30, phi = 40);title("truth")
par(old.par)
vis.scam(b,theta = 30, phi = 40)
## example with random effect smoother...
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)) # increasing smooth
f <- f1+f2
a <- factor(sample(1:10,200,replace=TRUE))
Xa <- model.matrix(~a-1) # random main effects
y <- f + Xa%*%rnorm(length(levels(a)))*.5 + rnorm(n)*.4
dat <- data.frame(x1=x1,x2=x2,y=y,a=a)
## fit model and plot...
b <- scam(y~s(x1,bs="cr")+s(x2,bs="mpi")+s(a,bs="re"), data=dat)
summary(b)
scam.check(b)
plot(b,pages=1,shade=TRUE)
## example with AR1 errors...
set.seed(8)
n <- 500
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)) # increasing smooth
f <- f1+f2
e <- rnorm(n,0,sd=2)
for (i in 2:n) e[i] <- .6*e[i-1] + e[i]
y <- f + e
dat <- data.frame(x1=x1,x2=x2,y=y)
b <- scam(y~s(x1,bs="cr")+s(x2,k=25,bs="mpi"),
data=dat, AR1.rho=.6, optimizer="efs")
b
## raw residuals still show correlation...
acf(residuals(b))
## but standardized are now fine...
acf(b$std.rsd)
}
Run the code above in your browser using DataLab