Merge branch 'jelena-markovic-logistic'

jonathan-taylor · jonathan-taylor · commit 69bf4732fedc · 2017-11-09T22:16:48.000-08:00
diff --git a/selectiveInference/R/funs.randomized.R b/selectiveInference/R/funs.randomized.R
@@ -6,6 +6,7 @@
 randomizedLasso = function(X, 
                            y, 
                            lam, 
+                           family=c("gaussian","binomial"),
                            noise_scale=NULL, 
                            ridge_term=NULL, 
                            noise_type=c('gaussian', 'laplace'),
@@ -17,6 +18,7 @@ randomizedLasso = function(X,
                            kkt_stop=TRUE,
                            parameter_stop=TRUE)
 {
+    family = match.arg(family)
 
     n = nrow(X); p = ncol(X)
     			
@@ -65,8 +67,8 @@ randomizedLasso = function(X,
     nactive = as.integer(0)
 
     result = solve_QP_wide(X,                  # design matrix
-    	                     lam / n,            # vector of Lagrange multipliers
-		                       ridge_term / n,     # ridge_term 
+                           lam / n,            # vector of Lagrange multipliers
+                           ridge_term / n,     # ridge_term 
                            max_iter, 
                            soln, 
                            linear_func, 
@@ -83,7 +85,7 @@ randomizedLasso = function(X,
 			                     parameter_stop)         # param_stop
     
     sign_soln = sign(result$soln)
-
+    
     unpenalized = lam == 0
     active = (!unpenalized) & (sign_soln != 0)
     inactive = (!unpenalized) & (sign_soln == 0)
@@ -110,8 +112,26 @@ randomizedLasso = function(X,
     I = inactive_set
     X_E = X[,E]
     X_I = X[,I]
-    L_E = t(X) %*% X[,E]
-
+    
+    if (length(E)==0){
+      return(list(active_set=c()))
+    }
+    
+    if (family=="binomial"){
+      unpen_reg = glm(y~X_E-1, family="binomial")
+      unpen_est = unpen_reg$coefficients
+      pi_fn = function(beta){
+        temp = X_E %*% as.matrix(beta)
+        return(as.vector(exp(temp)/(1+exp(temp)))) # n-dimensional
+      }
+      pi_vec = pi_fn(unpen_est)
+      W_E = diag(pi_vec*(1-pi_vec))
+    } else if (family=="gaussian"){
+      W_E = diag(rep(1,n))
+    }
+    
+    L_E = t(X) %*% W_E %*% X[,E]
+    
     coef_term = L_E
 
     signs_ = c(rep(1, sum(unpenalized)), sign_soln[active])
@@ -155,8 +175,12 @@ randomizedLasso = function(X,
                               offset_term = offset_term)
 
     # density for sampling optimization variables
-
+    
     observed_raw = -t(X) %*% y
+    if (family=="binomial"){
+      beta_E = result$soln[active_set]
+      observed_raw = observed_raw + t(X)%*%pi_fn(beta_E) - L_E %*% beta_E
+    }
     inactive_lam = lam[inactive_set]
     inactive_start = sum(unpenalized) + sum(active)
     active_start = sum(unpenalized)
@@ -191,11 +215,11 @@ randomizedLasso = function(X,
                 optimization_transform = opt_transform,
                 internal_transform = internal_transform,
                 log_optimization_density = log_optimization_density,
-		            observed_opt_state = observed_opt_state,
+                observed_opt_state = observed_opt_state,
                 observed_raw = observed_raw,
-		            noise_scale = noise_scale,
-		            soln = result$soln,
-		            perturb = perturb_
+                noise_scale = noise_scale,
+                soln = result$soln,
+                perturb = perturb_
                 ))
 
 }
@@ -330,6 +354,7 @@ conditional_density = function(noise_scale, lasso_soln) {
 randomizedLassoInf = function(X, 
                               y, 
                               lam, 
+                              family=c("gaussian", "binomial"),
                               sigma=NULL, 
                               noise_scale=NULL, 
                               ridge_term=NULL, 
@@ -349,10 +374,13 @@ randomizedLassoInf = function(X,
 
   n = nrow(X)
   p = ncol(X)
-
+  
+  family = match.arg(family)
+  
   lasso_soln = randomizedLasso(X, 
                                y, 
                                lam, 
+                               family=family,
                                noise_scale=noise_scale, 
                                ridge_term=ridge_term,
                                max_iter=max_iter,
@@ -409,17 +437,30 @@ randomizedLassoInf = function(X,
   X_E = X[, active_set]
   X_minusE = X[, inactive_set]
 
-  # if no sigma given, use OLS estimate
-
+  if (family == "gaussian") {
+    lm_y = lm(y ~ X_E - 1)
+    sigma_resid = sqrt(sum(resid(lm_y)^2) / lm_y$df.resid)
+    observed_target = lm_y$coefficients
+    W_E = diag(rep(1,n))
+    observed_internal = c(observed_target, t(X_minusE) %*% (y-X_E%*% observed_target))
+  } else if (family == "binomial") {
+    glm_y = glm(y~X_E-1)
+    sigma_resid = sqrt(sum(resid(glm_y)^2) / glm_y$df.resid)
+    observed_target = as.matrix(glm_y$coefficients)
+    temp = X_E%*%observed_target
+    pi_vec = exp(temp)/(1+exp(temp))
+    observed_internal =  c(observed_target, t(X_minusE) %*% (y-pi_vec))
+    W_E=diag(as.vector(pi_vec *(1-pi_vec)))
+  }
+  
+  # if no sigma given, use the estimate
+  
   if (is.null(sigma)) {
-        lm_y = lm(y ~ X_E - 1)
-        sigma = sqrt(sum(resid(lm_y)^2) / lm_y$df.resid)
+    sigma = sigma_resid
   }        
-
-  target_cov = solve(t(X_E) %*% X_E)*sigma^2
+  
+  target_cov = solve(t(X_E) %*% W_E %*% X_E)*sigma^2
   cov_target_internal = rbind(target_cov, matrix(0, nrow=p-nactive, ncol=nactive))
-  observed_target = solve(t(X_E) %*% X_E) %*% t(X_E) %*% y
-  observed_internal = c(observed_target, t(X_minusE) %*% (y-X_E%*% observed_target))
   internal_transform = lasso_soln$internal_transform
   opt_transform = lasso_soln$optimization_transform
   observed_raw = lasso_soln$observed_raw
@@ -494,5 +535,10 @@ randomizedLassoInf = function(X,
   return(list(active_set=active_set, pvalues=pvalues, ci=ci))
 }
 
+   
+    
+    
+    
+
 
 
diff --git a/selectiveInference/man/randomizedLassoInf.Rd b/selectiveInference/man/randomizedLassoInf.Rd
@@ -14,6 +14,7 @@ randomization.
 randomizedLassoInf(X, 
                    y, 
                    lam, 
+                   family=c("gaussian", "binomial"),
                    sigma=NULL, 
                    noise_scale=NULL, 
                    ridge_term=NULL, 
@@ -49,6 +50,9 @@ Value of lambda used to compute beta. See the above warning
  where obj is the object returned by glmnet (and [-1] removes the intercept,
  which glmnet always puts in the first component)
 } 
+\item{family}{
+Response type: "gaussian" (default), "binomial".
+}
 \item{sigma}{
 Estimate of error standard deviation. If NULL (default), this is estimated 
 using the mean squared residual of the full least squares based on 
diff --git a/tests/randomized/test_instances.R b/tests/randomized/test_instances.R
@@ -1,7 +1,7 @@
 library(selectiveInference)
 
-gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
-                             random_signs=TRUE, scale=TRUE, center=TRUE, seed=NA){
+get_instance = function(n, p, s, sigma=1, rho=0, signal=6, family="gaussian",
+                        X=NA, random_signs=TRUE, scale=TRUE, center=TRUE, seed=NA){
   if (!is.na(seed)){
     set.seed(seed)
   }
@@ -19,11 +19,20 @@ gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
     signs = sample(c(-1,1), s, replace = TRUE)
     beta = beta * signs
   }
-  y = X %*% beta + rnorm(n)*sigma
+  mu = X %*% beta
+  if (family=="gaussian"){
+    y = mu + rnorm(n)*sigma
+  } else if (family=="binomial"){
+    prob = exp(mu)/(1+exp(mu))
+    y= rbinom(n,1, prob)
+  }
   result = list(X=X,y=y,beta=beta)
   return(result)
 }
 
+
+
+
 test_randomized_lasso = function(n=100,p=200,s=0){
   set.seed(1)
   data = gaussian_instance(n=n,p=p,s=s, rho=0.3, sigma=3)
@@ -61,27 +70,29 @@ test_KKT=function(){
   
 
 
-collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=FALSE, lam=1.2){
+collect_results = function(n,p,s, nsim=100, level=0.9, 
+                           family = "binomial",
+                           condition_subgrad=FALSE, lam=1.2){
 
   rho=0.3
   sigma=1
   sample_pvalues = c()
   sample_coverage = c()
   for (i in 1:nsim){
-    data = gaussian_instance(n=n,p=p,s=s, rho=rho, sigma=sigma)
+    data = get_instance(n=n,p=p,s=s, rho=rho, sigma=sigma, family=family)
     X=data$X
     y=data$y
-    beta=data$beta
     result = selectiveInference:::randomizedLassoInf(X, y, 
-                                                     lam=lam, 
+                                                     lam, 
+                                                     family = family,
+                                                     sampler = "A",
                                                      sigma=sigma,
                                                      level=level, 
-                                                     sampler = "A",
                                                      burnin=1000, 
                                                      nsample=5000, 
                                                      condition_subgrad=condition_subgrad)
     if (length(result$active_set)>0){
-      true_beta = beta[result$active_set]
+      true_beta = data$beta[result$active_set]
       coverage = rep(0, nrow(result$ci))
       for (i in 1:nrow(result$ci)){
         if (result$ci[i,1]<true_beta[i] & result$ci[i,2]>true_beta[i]){
@@ -104,7 +115,7 @@ collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=FALSE,
 }
 
 set.seed(1)
-collect_results(n=100, p=20, s=0, lam=1.2)
+collect_results(n=100, p=20, s=0, lam=0.8)
 #test_randomized_lasso()
 #test_KKT()
 
