Nothing
## ---- echo = FALSE, results='hide', include=FALSE------------------------
library(mdpeer)
library(glmnet)
library(ggplot2)
## ---- fig.width = 3.3, fig.height = 2.8----------------------------------
library(mdpeer)
n <- 100
p1 <- 10
p2 <- 40
p <- p1 + p2
# Define graph adjacency matrix
A <- matrix(rep(0, p*p), nrow = p, ncol = p)
A[1:p1, 1:p1] <- 1
A[(p1+1):p, (p1+1):p] <- 1
diag(A) <- rep(0,p)
# Compute graph Laplacian matrix
L <- Adj2Lap(A)
vizu.mat(A, "graph adjacency matrix", 9); vizu.mat(L, "graph Laplacian matrix", 9)
## ------------------------------------------------------------------------
# simulate data objects
set.seed(1)
Z <- scale(matrix(rnorm(n*p), nrow = n, ncol = p))
X <- cbind(1, scale(matrix(rnorm(n*3), nrow = n, ncol = 3))) # cbind with 1s col. for intercept
b.true <- c(rep(1,p1), rep(0,p2)) # coefficients of interest (penalized)
beta.true <- c(0, runif(3)) # intercept=0 and 3 coefficients (non-penalized)
eta <- Z %*% b.true + X %*% beta.true
R2 <- 0.5 # assumed variance explained
sd.eps <- sqrt(var(eta) * (1 - R2) / R2)
error <- rnorm(n, sd = sd.eps)
y <- eta + error
## ------------------------------------------------------------------------
# estimate with riPEER
riPEER.out <- riPEER(Q = L, y = y, Z = Z, X = X, verbose = FALSE)
b.est.riPEER <- riPEER.out$b.est
# (graph penalty regulatization parameter, ridge penalty regularization parameter)
c(riPEER.out$lambda.Q, riPEER.out$lambda.R)
# estimate with cv.glmnet
library(glmnet)
cv.out.glmnet <- cv.glmnet(x = cbind(X,Z), y = y, alpha = 0, intercept = FALSE)
b.est.cv.glmnet <- unlist(coef(cv.out.glmnet))[5:54] # exclude intercept and covs
## ---- echo = FALSE, fig.width=7, fig.height=3----------------------------
plt.df <- data.frame(x = rep(1:p,2),
b.coef = c(b.est.riPEER,b.est.cv.glmnet),
est.method = c(rep("riPEER",p), rep("cv.glmnet",p)))
ggplot(plt.df, aes(x=x,y=b.coef,group=est.method, color=est.method)) +
geom_line(data = data.frame(x=1:p,y=b.true), aes(x=x,y=y),inherit.aes=FALSE) +
geom_line() + geom_point() + theme_bw(base_size = 9) +
labs(title = "black line: b.true", color = "estimation\nmethod")
## ------------------------------------------------------------------------
MSEr <- function(b.true, b.est) sum((b.true-b.est)^2)/sum(b.true^2)
# (riPEER error, cv.glmnet error)
c(MSEr(b.true,b.est.riPEER), MSEr(b.true,b.est.cv.glmnet))
## ------------------------------------------------------------------------
# simulate data objects
set.seed(1)
Z <- scale(matrix(rnorm(n*p), nrow = n, ncol = p))
X <- cbind(1, scale(matrix(rnorm(n*3), nrow = n, ncol = 3)))
b.true <- rep(c(-1,1), p/2)
# coefficients of interest (penalized)
beta.true <- c(0, runif(3)) # intercept=0 and 3 coefficients (non-penalized)
eta <- Z %*% b.true + X %*% beta.true
R2 <- 0.5 # assumed variance explained
sd.eps <- sqrt(var(eta) * (1 - R2) / R2)
error <- rnorm(n, sd = sd.eps)
y <- eta + error
## ------------------------------------------------------------------------
# estimate with riPEER
riPEER.out <- riPEER(Q = L, y = y, Z = Z, X = X, verbose = FALSE)
b.est.riPEER <- riPEER.out$b.est
c(riPEER.out$lambda.Q, riPEER.out$lambda.R)
# estimate with cv.glmnet
cv.out.glmnet <- cv.glmnet(x = cbind(X,Z), y = y, alpha = 0, intercept = FALSE)
b.est.cv.glmnet <- unlist(coef(cv.out.glmnet))[5:54] # exclude intercept and covs
## ---- echo = FALSE, fig.width=7, fig.height=3----------------------------
plt.df <- data.frame(x = rep(1:p,2),
b.coef = c(b.est.riPEER,b.est.cv.glmnet),
est.method = c(rep("riPEER",p), rep("cv.glmnet",p)))
ggplot(plt.df, aes(x=x,y=b.coef,group=est.method, color=est.method)) +
geom_line(data = data.frame(x=1:p,y=b.true), aes(x=x,y=y),inherit.aes=FALSE) +
geom_line() + geom_point() + theme_bw(base_size = 9) +
labs(title = "black line: b.true", color = "estimation\nmethod")
## ------------------------------------------------------------------------
MSEr <- function(b.true, b.est) sum((b.true-b.est)^2)/sum(b.true^2)
# (riPEER error, cv.glmnet error)
c(MSEr(b.true,b.est.riPEER), MSEr(b.true,b.est.cv.glmnet))
## ------------------------------------------------------------------------
# example based on `glmnet::cv.glmnet` CRAN documentation
set.seed(1010)
n=1000;p=100
nzc=trunc(p/10)
x=matrix(rnorm(n*p),n,p)
beta=rnorm(nzc)
fx= x[,seq(nzc)] %*% beta
eps=rnorm(n)*5
y=drop(fx+eps)
set.seed(1011)
cvob1=cv.glmnet(x,y, alpha = 0)
est.cv.glmnet <- coef(cvob1)[-c(1)] # exclude intercept and covs
# use riPEERc (X has column of 1s to represent intercept in a model)
riPEERc.out <- riPEERc(Q = diag(p), y = matrix(y), Z = x, X = matrix(rep(1,n)))
est.riPEERc <- riPEERc.out$b.est
## ---- echo = FALSE, fig.width=7, fig.height=3----------------------------
plt.df <- data.frame(x = rep(1:p,2),
b.coef = c(est.riPEERc,est.cv.glmnet),
est.method = c(rep("riPEERc",p), rep("cv.glmnet",p)))
ggplot(plt.df, aes(x=x,y=b.coef,group=est.method, color=est.method)) +
geom_line(data = data.frame(x=1:p,y=beta), aes(x=x,y=y),inherit.aes=FALSE) +
geom_line() + geom_point() + theme_bw(base_size = 9) +
labs(title = "black line: true", color = "estimation\nmethod")
## ------------------------------------------------------------------------
MSEr <- function(b.true, b.est) sum((b.true-b.est)^2)/sum(b.true^2)
# (riPEER error, cv.glmnet error)
c(MSEr(beta,est.riPEERc), MSEr(beta,est.cv.glmnet))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.