if (FALSE) {
## Concave P-splines example
## simulating data...
require(scam)
set.seed(1)
n <- 100
x <- sort(2*runif(n)-1)
f <- -4*x^2
y <- f + rnorm(n)*0.45
dat <- data.frame(x=x,y=y)
b <- scam(y~s(x,k=15,bs="cv"),family=gaussian,data=dat,not.exp=FALSE)
## fit unconstrained model...
b1 <- scam(y~s(x,k=15,bs="cr"),family=gaussian, data=dat,not.exp=FALSE)
## plot results ...
plot(x,y,xlab="x",ylab="y",cex=.5)
lines(x,f) ## the true function
lines(x,b$fitted,col=2) ## constrained fit
lines(x,b1$fitted,col=3) ## unconstrained fit
## Poisson version...
y <- rpois(n,15*exp(f))
dat <- data.frame(x=x,y=y)
## fit model ...
b <- scam(y~s(x,k=15,bs="cv"),family=poisson(link="log"),data=dat,not.exp=FALSE)
## fit unconstrained model...
b1<-scam(y~s(x,k=15,bs="cr"),family=poisson(link="log"), data=dat,not.exp=FALSE)
## plot results ...
plot(x,y,xlab="x",ylab="y",cex=.5)
lines(x,15*exp(f)) ## the true function
lines(x,b$fitted,col=2) ## constrained fit
lines(x,b1$fitted,col=3) ## unconstrained fit
## plotting on log scale...
plot(x,log(15*exp(f)),type="l",cex=.5) ## the true function
lines(x,log(b$fitted),col=2) ## constrained fit
lines(x,log(b1$fitted),col=3) ## unconstrained fit
## 'by' factor example...
set.seed(9)
n <- 400
x <- sort(runif(n,-.5,.5))
f1 <- -.7*x+cos(x)-3
f2 <- -20*x^2
par(mfrow=c(1,2))
plot(x,f1,type="l");plot(x,f2,type="l")
e <- rnorm(n, 0, 1.5)
fac <- as.factor(sample(1:2,n,replace=TRUE))
fac.1 <- as.numeric(fac==1)
fac.2 <- as.numeric(fac==2)
y <- f1*fac.1 + f2*fac.2 + e
dat <- data.frame(y=y,x=x,fac=fac,f1=f1,f2=f2)
b2 <- scam(y ~ fac+s(x,by=fac,bs="cv"),data=dat,optimizer="efs")
plot(b2,pages=1,scale=0,shade=TRUE)
summary(b2)
x11()
vis.scam(b2,theta=50,color="terrain")
## numeric 'by' variable example...
set.seed(6)
n <- 100
x <- sort(2*runif(n)-1)
z <- runif(n,-2,3)
f <- -4*x^2
y <- f*z + rnorm(n)*0.6
dat <- data.frame(x=x,z=z,y=y)
b <- scam(y~s(x,k=15,by=z,bs="cvBy"),data=dat)
summary(b)
par(mfrow=c(1,2))
plot(b,shade=TRUE)
## unconstrained fit...
b1 <- scam(y~s(x,k=15,by=z),data=dat)
plot(b1,shade=TRUE)
summary(b1)
}
Run the code above in your browser using DataLab