library(survival)

set.seed(23)
n <- 250

X <- runif(n, 1, 11) # uniform on (1,11)
EX <- log(6)
Z <- rnorm(n, 0, 1) # standard Gaussian
Y <- log(X) - log(6) + Z

q <- c(1,3,5,7,9,11)
qGroups <- cut(X, q)
qMean <- tapply(X,qGroups,mean)

x0 <- ifelse(X <= 3, 1,0)
x1 <- ifelse(3< X & X <= 5, 1, 0)
x2 <- ifelse(5< X & X <= 7, 1, 0)
x3 <- ifelse(7< X & X <= 9, 1, 0)
x4 <- ifelse(9 < X, 1,0)
xstar <- x0*2+x1*4+x2*6+x3*8+x4*10
x1234 <- x1+x2+x3+x4
x234 <- x2+x3+x4
x34 <- x3+x4
x1plus <- x1234*(X-3)
x2plus <- x234*(X-5)
x3plus <- x34*(X-7)
x4plus <- x4*(X-9)
xsq <- X*X;
x1plussq <- x1plus*x1plus
x2plussq <- x2plus*x2plus
x3plussq <- x3plus*x3plus
x4plussq <- x4plus*x4plus

x <- matrix(nrow=length(X), ncol=5)
qx <- matrix(nrow=length(X), ncol=5)
qrx <- matrix(nrow=length(X), ncol=4)

model1 <- lm(Y ~ x0 + x1 + x2 + x3 + x4 - 1)
model2 <- lm(Y ~ X + x1plus + x2plus + x3plus + x4plus)
model3 <- lm(Y ~ X + xsq + x1plussq + x2plussq + x3plussq + x4plussq)

## The functions are 'vectorized', so that they can take an vector as argument
## This is needed for use in curve()

linpredSpline <- function(x){
  model <- model2
  f <- function(x){
    y <- model$coefficients[1] + model$coefficients[2]*(x)
    for (i in c(3:6)){
      y <- y + model$coefficients[i]*ifelse(x > q[i-1], x-q[i-1], 0)
    }
    return(y)
  }
  ## here the internal function, f, is applied to each entry in x.
  sapply(x, f)
}

newX <- seq(1, 11, 0.01)

linpredQuadratic <- function(x){
  model <- model3
  f <- function(x){
    y <- model$coefficients[1] + model$coefficients[2]*(x) +
      model$coefficients[3]*(x^2) 
    for (i in c(4:7)){
      y <- y + model$coefficients[i]*ifelse(x > q[i-2], (x-q[i-2])^2, 0)
    }
    return(y) 
  }
  sapply(x, f)
}

basecurve <- function(...){
  curve(log(x)-log(6),
        from = 1,
        to = 11,
        n = 2000,
        lty = "dotted"
        ,...)
}

par(mfrow=c(2,2))
basecurve(
          ylim = c(-2,.75),
          xlab = "x",
          ylab = "Linear Predictor")
lines(predict(lm(Y~X), data.frame(X = newX)) ~ newX)

plot(q, c(model1$coefficients, model1$coefficients[5]) ,
     type="s",
     xlim = c(1,11),
     ylim = c(-2,.75),
     xlab = "x",
     ylab = "Linear Predictor")

basecurve(add=TRUE)

curve(linpredSpline(x),
      from = 1,
      to = 11,
      ##        yaxp = c(-1,5,6),
      ylim = c(-2,.75),
      xlab = "x",
      ylab = "Linear predictor",
      lty = "solid")
basecurve(add=TRUE)

curve(linpredQuadratic(x),
      from = 1,
      to = 11,
      ylim = c(-2,.75),
      xlab = "x",
      ylab = "Linear Predictor",
      n=2000,
      lty = "solid")
basecurve(add=TRUE)