# requires "care" R package version 1.1.1

library("care")

# load Lu et al. (2004) data set
data(lu2004)
x = lu2004$x
y = lu2004$y # age
dim(x) # 30 403

# regularization parameters and number of included predictors
reg = slm(x,y)$regularization
lambda = reg[1] # 0.1373293
lambda.var = reg[2] # 0.0245981
numpred =  c(seq(5, 45, by=10), seq(53, 403, by=25))
numpred

# ordering according to CAR scores 
car = carscore(x, y, lambda=lambda)
ocar = order(car^2, decreasing=TRUE)
car.predlist = make.predlist(ocar, numpred, name="CAR")
car.models = slm.models(x, y, car.predlist, lambda=lambda, lambda.var=lambda.var)

# ordering accoding to marginal correlations 
marg = carscore(x, y, lambda=lambda, diagonal=TRUE) # shrinking not actually reqired
omarg = order(marg^2, decreasing=TRUE)
marg.predlist = make.predlist(omarg, numpred, name="MARG")
marg.models = slm.models(x, y, marg.predlist, lambda=lambda, lambda.var=lambda.var)

# make plot
ylim = range( c(marg.models$R2, car.models$R2) )
plot(car.models$numpred, car.models$R2, type="b", ylim=ylim,
  xlab="number of included predictors",
  ylab="model fit (R2)",
  main="CAR and Marginal Correlation Models for Brain Data")
points(marg.models$numpred, marg.models$R2, col=2, type="b")
legend("bottomright", c("CAR Score", "Marginal Correlation"), col=c(1,2), lty=c(1,1) )