toc

Author: zuoyi93

.Rhistory

@@ -0,0 +1,86 @@
+usethis::use_vignette("my-vignette")
+citation("brglm2")
+library(devtools)
+library(roxygen2)
+document()
+check()
+x=1:10
+y=x*2+1+rnorm(10)
+lm(y~x)
+lm.fit(y~x)
+.lm.fit(y~x)
+.lm.fit(x=x,y=y)
+cat("hi",y,"\n")
+document()
+check()
+install()
+document()
+check()
+document()
+check()
+document()
+check()
+sign(1)
+sign(2)
+sign(-2)
+document()
+document()
+check()
+load_all()
+data.linear <- gen.sim.data()
+abs(0)
+document()
+check()
+data.linear <- gen.sim.data()
+head(data.linear)
+0^2
+0^0
+rho=-0.2
+p=4
+cov.structure <- matrix(0, p, p)
+for (i in 1:p) {
+for (j in 1:p) {
+if(rho != 0){
+cov.structure[i, j] <- sign(rho) * abs(rho)^(abs(i - j))
+}else{
+cov.structure[i, j] <- abs(rho)^(abs(i - j))
+}
+}
+}
+cov.structure
+if(rho < 0) diag(cov.structure) <- 1
+cov.structure
+data.linear <- gen.sim.data(rho=-0.2,n=20,p=4,s=2)
+load_all()
+data.linear <- gen.sim.data(rho=-0.2,n=20,p=4,s=2)
+data.linear[[1]]
+cov(data.linear[[1]])
+cor(data.linear[[1]])
+data.linear <- gen.sim.data(rho=-0.2,n=40,p=4,s=2)
+cor(data.linear[[1]])
+data.linear <- gen.sim.data(rho=-0.5,n=80,p=4,s=2)
+data.linear <- gen.sim.data(rho=-0.4,n=80,p=4,s=2)
+cor(data.linear[[1]])
+document()
+check()
+data.linear <- gen.sim.data(n = 20)
+x <- data.linear[[1]]
+head(x)
+data.linear <- gen.sim.data(n = 20, p = 10, s = 4)
+x <- data.linear[[1]]
+x
+data.linear <- gen.sim.data(n = 20, p = 10, s = 4, family="cox")
+data.linear[[2]]
+glmnet
+x
+y
+data.linear <- gen.sim.data(n = 20, p = 10, s = 4)
+x <- data.linear[[1]]
+y <- data.linear[[2]]
+index <- data.linear[[3]]
+true.beta <- data.linear[[4]]
+temp=glmnet(x,y)
+names(temp)
+document()
+document()
+getwd()

NAMESPACE

@@ -5,8 +5,9 @@ S3method(plot,sgpv)
 S3method(predict,sgpv)
 S3method(print,sgpv)
 S3method(summary,sgpv)
-export(gen.data)
+export(gen.sim.data)
 export(pro.sgpv)
+export(which.sgpv)
 importFrom(MASS,mvrnorm)
 importFrom(glmnet,cv.glmnet)
 importFrom(glmnet,glmnet)
@@ -19,4 +20,8 @@ importFrom(stats,coef)
 importFrom(stats,complete.cases)
 importFrom(stats,lm)
 importFrom(stats,predict)
+importFrom(stats,rbinom)
+importFrom(stats,rexp)
 importFrom(stats,rnorm)
+importFrom(stats,rpois)
+importFrom(stats,runif)

R/data.R

@@ -39,3 +39,33 @@
 #' }
 
 "t.housing"
+
+
+#' Spine data
+#'
+#' Lower back pain can be caused by a variety of problems with any parts of the complex,
+#' interconnected network of spinal muscles, nerves, bones, discs or tendons
+#' in the lumbar spine. This dataset contains 12 biomechanical attributes from
+#' 310 patients, of whom 100 are normal and 210 are abnormal (Disk Hernia or
+#' Spondylolisthesis). The goal is to differentiate the normal patients from the
+#' abnormal using those 12 variables.
+#'
+#' @source \url{http://archive.ics.uci.edu/ml/datasets/vertebral+column}
+#' @format
+#' \describe{
+#'  \item{pelvic_incidence}{pelvic incidence}
+#'  \item{pelvic_tilt}{pelvic tilt}
+#'  \item{lumbar_lordosis_angle}{lumbar lordosis angle}
+#'  \item{sacral_slope}{sacral slope}
+#'  \item{pelvic_radius}{pelvic radius}
+#'  \item{degree_spondylolisthesis}{degree of spondylolisthesis}
+#'  \item{pelvic_slope}{pelvic slope}
+#'  \item{direct_tilt}{direct tilt}
+#'  \item{thoracic_slope}{thoracic slope }
+#'  \item{cervical_tilt}{cervical tilt}
+#'  \item{sacrum_angle}{sacrum angle}
+#'  \item{scoliosis_slope}{scoliosis slope}
+#'  \item{outcome}{1 is abnormal (Disk Hernia or Spondylolisthesis) and 0 is normal}
+#' }
+
+"spine"

R/pro.sgpv.R

@@ -8,7 +8,11 @@
 #' @param x Independent variables, can be a \code{matrix} or a \code{data.frame}
 #' @param y Dependent variable, can be a \code{vector} or a column from a \code{data.frame}
 #' @param stage Algorithm indicator. 1 denotes the one-stage algorithm and
-#' 2 denotes the two-stage algorithm. Default is 1.
+#' 2 denotes the two-stage algorithm. Default is 2. When \code{n} is less than \code{p},
+#' only the two-stage algorithm is available.
+#' @param family A description of the error distribution and link function to be
+#'  used in the model. It can take the value of `\code{gaussian}`, `\code{binomial}`,
+#'  `\code{poisson}`, and `\code{cox}`. Default is `\code{gaussian}`
 #'
 #' @return A list of following components:
 #' \describe{
@@ -17,6 +21,8 @@
 #' \item{lambda}{Cross-validated lambda in the two-stage algorithm. \code{NULL} for the one-stage algorithm}
 #' \item{x}{Input data \code{x}}
 #' \item{y}{Input data \code{y}}
+#' \item{family}{\code{family} from the input}
+#' \item{stage}{\code{stage} from the input}
 #' }
 #' @export
 #' @seealso
@@ -27,19 +33,16 @@
 #' * [plot.sgpv()] plots variable selection results
 #' @examples
 #'
-#' # load the package
-#' library(ProSGPV)
-#'
 #' # prepare the data
 #' x <- t.housing[, -ncol(t.housing)]
 #' y <- t.housing$V9
 #'
-#' # run one-stage algorithm
-#' out.sgpv.1 <- pro.sgpv(x = x, y = y, stage = 1)
+#' # run ProSGPV in linear regression
+#' out.sgpv <- pro.sgpv(x = x, y = y)
 #'
-#' # More examples at https://github.com/zuoyi93/ProSGPV
-pro.sgpv <- function(x, y, stage = c(1, 2)) {
-  if (!(stage %in% 1:2)) stop("Stage only takes value of 1 or 2.")
+#' # More examples at https://github.com/zuoyi93/ProSGPV/tree/master/vignettes
+pro.sgpv <- function(x, y, stage = c(1, 2),
+                     family = c("gaussian", "binomial", "poisson", "cox")) {
 
   if (nrow(x) != length(y)) stop("Input x and y have different number of observations")
 
@@ -48,19 +51,46 @@ pro.sgpv <- function(x, y, stage = c(1, 2)) {
   if (any(complete.cases(x) == F) | any(complete.cases(y) == F)) {
     warning("Only complete records will be used.")
     comp.index <- complete.cases(data.frame(x, y))
-    x <- x[comp.index, ]
-    y <- y[comp.index]
+    if(family != "cox"){
+      x <- x[comp.index, ]
+      y <- y[comp.index]
+    }else{
+      x <- x[comp.index, ]
+      y <- y[comp.index,]
+    }
+
   }
 
+  if(missing(stage)) stage <- 2
+  if(missing(family)) family <- "gaussian"
+
+  stage <- match.arg(stage)
+  family <- match.arg(family)
+
+  if(stage == 1 & nrow(x)<ncol(x)) stage <- 2
+
   if (is.null(colnames(x))) colnames(x) <- paste("V", 1:ncol(x), sep = "")
 
-  xs <- scale(x)
-  ys <- scale(y)
+  if(family == "gaussian"){
+    xs <- scale(x)
+    ys <- scale(y)
+  }else{
+    xs <- as.matrix(x)
+    ys <- y
+  }
+
 
   if (stage == 2) {
-    lasso.cv <- cv.glmnet(xs, ys)
-    lambda <- lasso.cv$lambda.1se
-    candidate.index <- which(coef(lasso.cv, s = lambda)[-1] != 0)
+
+    if(family != "cox"){
+
+      lasso.cv <- cv.glmnet(xs, ys, family = family)
+      lambda <- lasso.cv$lambda.1se
+      candidate.index <- which(coef(lasso.cv, s = lambda)[-1] != 0)
+    }else{
+      candidate.index <- which(coef(lasso.cv, s = lambda) != 0)
+    }
+
   } else {
     candidate.index <- 1:ncol(xs)
     lambda <- NULL
@@ -72,8 +102,10 @@ pro.sgpv <- function(x, y, stage = c(1, 2)) {
     var.index = out.sgpv,
     var.label = colnames(x)[out.sgpv],
     lambda = lambda,
-    x = x,
-    y = y
+    x = data.frame(x),
+    y = y,
+    family = family,
+    stage = stage
   )
 
   class(out) <- "sgpv"
@@ -104,9 +136,9 @@ pro.sgpv <- function(x, y, stage = c(1, 2)) {
 #' out.sgpv.1
 print.sgpv <- function(x, ...) {
   if (length(x$var.index) > 0) {
-    cat("Selected variables are", x$var.label)
+    cat("Selected variables are", x$var.label, "\n")
   } else {
-    cat("None of variables are selected.")
+    cat("None of variables are selected.\n")
   }
 }