#define function that returns positive part
pos <- function(x) {
  res <- x
  res[x<0] <- 0
  res
}

n <- 1000
set.seed(63412)
x <- 2*rnorm(n)

#define cubic spline basis functions
a <- -1
b <- 0
c <- 1
x1 <- x
x2 <- x^2
x3 <- x^3
#knots at -1, 0, 1
x4 <- pos(x-(-1))^3
x5 <- pos(x-0)^3
x6 <- pos(x-1)^3

#specify true coefficients
beta <- c(0,1,0.1,0.2,-0.1,0.4,0.3)

y <- beta[1]+beta[2]*x1+beta[3]*x2+beta[4]*x3+beta[5]*x4+beta[6]*x5+beta[7]*x6+rnorm(n)

plot(x,y)

#make some x values missing MCAR
x[runif(n)<0.5] <- NA
newdata <- data.frame(x,y)

#impute using smcfcs
library(smcfcs)
imps <- smcfcs(newdata, smtype="lm", 
               smformula="y~x+I(x^2)+I(x^3)+I(pos(x-(-1))^3)+I(pos(x-0)^3)+I(pos(x-1)^3)",
               m=1, method=c("norm",""))
with(imps$impDatasets[[1]], plot(x,y))

#reimpute with larger rjlimit
imps <- smcfcs(newdata, smtype="lm", 
               smformula="y~x+I(x^2)+I(x^3)+I(pos(x-(-1))^3)+I(pos(x-0)^3)+I(pos(x-1)^3)",
               m=1, method=c("norm",""),rjlimit=100000)
with(imps$impDatasets[[1]], plot(x,y))