##=============================
## 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)
par(mfrow=c(2,2))
plot(b,shade=TRUE)
if (FALSE) {
## bivariate example...
## simulating data...
set.seed(2)
n <- 30
x1 <- sort(runif(n)*4-1)
x2 <- sort(runif(n))
f1 <- matrix(0,n,n)
for (i in 1:n) for (j in 1:n)
{ f1[i,j] <- -exp(4*x1[i])/(1+exp(4*x1[i]))+2*sin(pi*x2[j])}
f <- as.vector(t(f1))
y <- f+rnorm(length(f))*.2
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)
## fit model and plot ...
b <- scam(y~s(x1,x2,k=c(10,10),bs=c("tesmd1","ps")),data=dat,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)
plot(y,b$fitted.values,xlab="Simulated data",ylab="Fitted data",pch=".",cex=3)
par(old.par)
## 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...
x11()
acf(b$std.rsd)
}
Run the code above in your browser using DataLab