Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ???
To elaborate on Frank's response, the analysis plan of 1. Look at the data and select important variables 2. Put that truncated list into your favorite statistic procedure 3. Ask - are the p-values (c-statistic, coefficients, .) reliable? is a very old plan. The answer to the last question is always NO. The traditional step 1 involved graphs and t-tests, but replacing it by something modern does not change the overall outcome. Terry T __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ???
On Sep 29, 2013, at 2:16 PM, E Joffe wrote: HI, Thank you for your answer. There were 301 events out of 394 observations. Study goals: Identify proteins with prognostic power in patients with AML. There were 232 proteins studied. Traditional models won't converge. I wanted to do a multivariate survival analysis that would allow me to identify prognostic proteins and get an estimation of their hazard ratio. LASSO allowed me to model the data but in order to get the HRs I needed to run it a second time throw coxph - this was done based on a discussion on stackexchange http://r.789695.n4.nabble.com/estimating-survival-times-with-glmnet-and-coxph-td4614225.html which I tweaked a bit. The features selected by LASSO are not necessarily all significant when constructing a coxph model in particular due to the lack of penalization. This was commented on by Dr. Hastie in one of his papers (or in the video lecture he gave - I need to find the exact source). I understand that by doing this procedure I might be reducing the accuracy of the model which is why I used a bootstrap evaluation. Interestingly, the coxph models following stepwise feature selection had marginally better performance in the evaluation compared with the LASSO models (maybe the LASSO over-fit the data a bit ? - don't know) What I still don't understand is why with a C-Index of 0.76 and a Brier score of 0.07 when I look at the plot of the coxph object the confidence intervals are so far ? The item you referenced was not on stackexchange but rather on Nabble. The archived copy in R-help is here: [R] Estimating survival times with glmnet and coxph https://stat.ethz.ch/pipermail/r-help/2012-May/312029.html https://stat.ethz.ch/pipermail/r-help/2012-May/312038.html The response from Terry Therneau at that time explained why the CI's were wide. I think this implies that even without the additional (faulty) step of stepwise reduction in covariates you still would not be able to use the CI's coxphfrom such a fit. -- David. Further, I am wondering whether there is a sanity check I could perform to make sure there is no error in the Brier/C-Index evaluation. Thanks again for all you advice on this project, Erel -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Saturday, September 28, 2013 4:41 PM To: E Joffe Cc: r-help@r-project.org Subject: Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ??? On Sep 28, 2013, at 2:39 AM, E Joffe wrote: Hi all, I am using COX LASSO (glmnet / coxnet) regression to analyze a dataset of 394 obs. / 268 vars. I use the following procedure: 1. Construct a coxnet on the entire dataset (by cv.glmnet) 2. Pick the significant features by selecting the non-zero coefficient under the best lambda selected by the model 3. Build a coxph model with bi-directional stepwise feature selection limited to the coxnet selected features. I was a bit puzzled by the third step. Once you had a reduced model from glmnet, what was the statistical basis for further elimination of variables? (Quite apart from the statistical issues, I was rather surprised that this procedure even produced results since the 'step' function is not described in the 'stats' package as applying to 'coxph' model objects.) To validate the model I use both Brier score (library=peperr) and Harrel's [Harrell] C-Index (library=Hmisc) with a bootstrap of 50 iterations. MY QUESTION : I am getting fairly good C-Index (0.76) and Brier (0.07) values for the models however per the coxnet the %Dev explained by the model is at best 0.27 and when I plot the survfit of the coxph the plotted confidence interval is very large. How many events did you have? (The width of CI's is most importantly dependent on event counts and not particularly improved by a high case count. The power considerations are very similar to those of a binomial test.) What am I missing here ? Perhaps sufficient events? (You also seem to be missing a description of the study goals.) -- David. %DEV=27% Brier score - 0.07 ($ipec.coxglmnet - [1] 7.24) C-Index - 0.76 ($cIndex - [1] 0.763) DATA: [Private Health Information - can't publish] 394 obs./268 vars. CODE (need to define a dataset with 'time' and 'status' variables): library(survival) library (glmnet) library (c060) library (peperr) library (Hmisc) #creat Y (survival matrix) for glmnet surv_obj - Surv(dataset$time,dataset$status) ## tranform categorical variables into binary variables with dummy for dataset predict_matrix - model.matrix(~ ., data=dataset, contrasts.arg = lapply (dataset[,sapply(dataset, is.factor)], contrasts)) ## remove the statu/time variables from the predictor matrix (x) for glmnet
Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ???
HI, Thank you for your answer. There were 301 events out of 394 observations. Study goals: Identify proteins with prognostic power in patients with AML. There were 232 proteins studied. Traditional models won't converge. I wanted to do a multivariate survival analysis that would allow me to identify prognostic proteins and get an estimation of their hazard ratio. LASSO allowed me to model the data but in order to get the HRs I needed to run it a second time throw coxph - this was done based on a discussion on stackexchange http://r.789695.n4.nabble.com/estimating-survival-times-with-glmnet-and-coxp h-td4614225.html which I tweaked a bit. The features selected by LASSO are not necessarily all significant when constructing a coxph model in particular due to the lack of penalization. This was commented on by Dr. Hastie in one of his papers (or in the video lecture he gave - I need to find the exact source). I understand that by doing this procedure I might be reducing the accuracy of the model which is why I used a bootstrap evaluation. Interestingly, the coxph models following stepwise feature selection had marginally better performance in the evaluation compared with the LASSO models (maybe the LASSO over-fit the data a bit ? - don't know) What I still don't understand is why with a C-Index of 0.76 and a Brier score of 0.07 when I look at the plot of the coxph object the confidence intervals are so far ? Further, I am wondering whether there is a sanity check I could perform to make sure there is no error in the Brier/C-Index evaluation. Thanks again for all you advice on this project, Erel -Original Message- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Saturday, September 28, 2013 4:41 PM To: E Joffe Cc: r-help@r-project.org Subject: Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ??? On Sep 28, 2013, at 2:39 AM, E Joffe wrote: Hi all, I am using COX LASSO (glmnet / coxnet) regression to analyze a dataset of 394 obs. / 268 vars. I use the following procedure: 1.Construct a coxnet on the entire dataset (by cv.glmnet) 2.Pick the significant features by selecting the non-zero coefficient under the best lambda selected by the model 3.Build a coxph model with bi-directional stepwise feature selection limited to the coxnet selected features. I was a bit puzzled by the third step. Once you had a reduced model from glmnet, what was the statistical basis for further elimination of variables? (Quite apart from the statistical issues, I was rather surprised that this procedure even produced results since the 'step' function is not described in the 'stats' package as applying to 'coxph' model objects.) To validate the model I use both Brier score (library=peperr) and Harrel's [Harrell] C-Index (library=Hmisc) with a bootstrap of 50 iterations. MY QUESTION : I am getting fairly good C-Index (0.76) and Brier (0.07) values for the models however per the coxnet the %Dev explained by the model is at best 0.27 and when I plot the survfit of the coxph the plotted confidence interval is very large. How many events did you have? (The width of CI's is most importantly dependent on event counts and not particularly improved by a high case count. The power considerations are very similar to those of a binomial test.) What am I missing here ? Perhaps sufficient events? (You also seem to be missing a description of the study goals.) -- David. %DEV=27% Brier score - 0.07 ($ipec.coxglmnet - [1] 7.24) C-Index - 0.76 ($cIndex - [1] 0.763) DATA: [Private Health Information - can't publish] 394 obs./268 vars. CODE (need to define a dataset with 'time' and 'status' variables): library(survival) library (glmnet) library (c060) library (peperr) library (Hmisc) #creat Y (survival matrix) for glmnet surv_obj - Surv(dataset$time,dataset$status) ## tranform categorical variables into binary variables with dummy for dataset predict_matrix - model.matrix(~ ., data=dataset, contrasts.arg = lapply (dataset[,sapply(dataset, is.factor)], contrasts)) ## remove the statu/time variables from the predictor matrix (x) for glmnet predict_matrix - subset (predict_matrix, select=c(-time,-status)) ## create a glmnet cox object using lasso regularization and cross validation glmnet.cv - cv.glmnet (predict_matrix, surv_obj, family=cox) ## get the glmnet model on the full dataset glmnet.obj - glmnet.cv$glmnet.fit # find lambda index for the models with least partial likelihood deviance (by cv.glmnet) optimal.lambda - glmnet.cv$lambda.min# For a more parsimoneous model use lambda.1se lambda.index - which(glmnet.obj$lambda==optimal.lambda) # take beta for optimal lambda optimal.beta - glmnet.obj$beta[,lambda.index] # find non zero beta coef
Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ???
This appears to be a statistics, not an R-help question, so should
probably be asked on a statistics list, not here (e.g.
stats.stackexchange.com).
But if I understand your issue correctly, perhaps the heart f the
matter is: why do you think a stable fit must explain a lot of the
variation? Feel free to ignore if I I'm wrong or discuss further on
the statistics list.
Cheers,
Bert
On Sat, Sep 28, 2013 at 2:39 AM, E Joffe ejo...@hotmail.com wrote:
Hi all,
I am using COX LASSO (glmnet / coxnet) regression to analyze a dataset of
394 obs. / 268 vars.
I use the following procedure:
1. Construct a coxnet on the entire dataset (by cv.glmnet)
2. Pick the significant features by selecting the non-zero coefficient
under the best lambda selected by the model
3. Build a coxph model with bi-directional stepwise feature selection
limited to the coxnet selected features.
To validate the model I use both Brier score (library=peperr) and Harrel's
C-Index (library=Hmisc) with a bootstrap of 50 iterations.
MY QUESTION : I am getting fairly good C-Index (0.76) and Brier (0.07)
values for the models however per the coxnet the %Dev explained by the model
is at best 0.27 and when I plot the survfit of the coxph the plotted
confidence interval is very large.
What am I missing here ?
%DEV=27%
Brier score - 0.07 ($ipec.coxglmnet - [1] 7.24)
C-Index - 0.76 ($cIndex - [1] 0.763)
DATA: [Private Health Information - can't publish] 394 obs./268 vars.
CODE (need to define a dataset with 'time' and 'status' variables):
library(survival)
library (glmnet)
library (c060)
library (peperr)
library (Hmisc)
#creat Y (survival matrix) for glmnet
surv_obj - Surv(dataset$time,dataset$status)
## tranform categorical variables into binary variables with dummy for
dataset
predict_matrix - model.matrix(~ ., data=dataset,
contrasts.arg = lapply
(dataset[,sapply(dataset, is.factor)], contrasts))
## remove the statu/time variables from the predictor matrix (x) for
glmnet
predict_matrix - subset (predict_matrix, select=c(-time,-status))
## create a glmnet cox object using lasso regularization and cross
validation
glmnet.cv - cv.glmnet (predict_matrix, surv_obj, family=cox)
## get the glmnet model on the full dataset
glmnet.obj - glmnet.cv$glmnet.fit
# find lambda index for the models with least partial likelihood
deviance (by cv.glmnet)
optimal.lambda - glmnet.cv$lambda.min# For a more parsimoneous
model use lambda.1se
lambda.index - which(glmnet.obj$lambda==optimal.lambda)
# take beta for optimal lambda
optimal.beta - glmnet.obj$beta[,lambda.index]
# find non zero beta coef
nonzero.coef - abs(optimal.beta)0
selectedBeta - optimal.beta[nonzero.coef]
# take only covariates for which beta is not zero
selectedVar - predict_matrix[,nonzero.coef]
# create a dataframe for trainSet with time, status and selected
variables in binary representation for evaluation in pec
reformat_dataSet - as.data.frame(cbind(surv_obj,selectedVar))
# glmnet.cox only with meaningful features selected by stepwise
bidirectional AIC feature selection
glmnet.cox.meaningful - step(coxph(Surv(time,status) ~
.,data=reformat_dataSet),direction=both)
##--
-
##MODEL PERFORMANCE
##--
-
##
## Calculate the Brier score - pec does its own bootstrap so this
function runs on i=51 (i.e., whole trainset)
## Brier score calculation to cox-glmnet
peperr.coxglmnet - peperr(response=surv_obj, x=selectedVarCox,
fit.fun=fit.coxph, load.all=TRUE,
indices=resample.indices(n=nrow(surv_obj),
method=boot, sample.n=50))
# Get error predictions from peperr
prederr.coxglmnet - perr(peperr.coxglmnet)
# Integrated prediction error Brier score calculation
ipec.coxglmnet-ipec(prederr.coxglmnet,
eval.times=peperr.coxglmnet$attribute, response=surv_obj)
## C-Index calculation 50 iter bootstrapping
for (i in 1:50){
print (paste(Iteration:,i))
train - sample(1:nrow(dataset), nrow(dataset), replace = TRUE) ##
random sampling with replacement
# create a dataframe for trainSet with time, status and selected
variables in binary representation for evaluation in pec
reformat_trainSet - reformat_dataSet [train,]
# glmnet.cox only with meaningful features selected by stepwise
bidirectional AIC feature selection
glmnet.cox.meaningful.test - step(coxph(Surv(time,status) ~
.,data=reformat_dataSet),direction=both)
selectedVarCox -
Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ???
On Sep 28, 2013, at 2:39 AM, E Joffe wrote:
Hi all,
I am using COX LASSO (glmnet / coxnet) regression to analyze a
dataset of
394 obs. / 268 vars.
I use the following procedure:
1. Construct a coxnet on the entire dataset (by cv.glmnet)
2. Pick the significant features by selecting the non-zero coefficient
under the best lambda selected by the model
3. Build a coxph model with bi-directional stepwise feature selection
limited to the coxnet selected features.
I was a bit puzzled by the third step. Once you had a reduced model
from glmnet, what was the statistical basis for further elimination of
variables?
(Quite apart from the statistical issues, I was rather surprised that
this procedure even produced results since the 'step' function is not
described in the 'stats' package as applying to 'coxph' model objects.)
To validate the model I use both Brier score (library=peperr) and
Harrel's [Harrell]
C-Index (library=Hmisc) with a bootstrap of 50 iterations.
MY QUESTION : I am getting fairly good C-Index (0.76) and Brier
(0.07)
values for the models however per the coxnet the %Dev explained by
the model
is at best 0.27 and when I plot the survfit of the coxph the plotted
confidence interval is very large.
How many events did you have? (The width of CI's is most importantly
dependent on event counts and not particularly improved by a high case
count. The power considerations are very similar to those of a
binomial test.)
What am I missing here ?
Perhaps sufficient events? (You also seem to be missing a description
of the study goals.)
--
David.
%DEV=27%
Brier score - 0.07 ($ipec.coxglmnet - [1] 7.24)
C-Index - 0.76 ($cIndex - [1] 0.763)
DATA: [Private Health Information - can't publish] 394 obs./268 vars.
CODE (need to define a dataset with 'time' and 'status' variables):
library(survival)
library (glmnet)
library (c060)
library (peperr)
library (Hmisc)
#creat Y (survival matrix) for glmnet
surv_obj - Surv(dataset$time,dataset$status)
## tranform categorical variables into binary variables with
dummy for
dataset
predict_matrix - model.matrix(~ ., data=dataset,
contrasts.arg = lapply
(dataset[,sapply(dataset, is.factor)], contrasts))
## remove the statu/time variables from the predictor matrix (x)
for
glmnet
predict_matrix - subset (predict_matrix, select=c(-time,-status))
## create a glmnet cox object using lasso regularization and cross
validation
glmnet.cv - cv.glmnet (predict_matrix, surv_obj, family=cox)
## get the glmnet model on the full dataset
glmnet.obj - glmnet.cv$glmnet.fit
# find lambda index for the models with least partial likelihood
deviance (by cv.glmnet)
optimal.lambda - glmnet.cv$lambda.min# For a more parsimoneous
model use lambda.1se
lambda.index - which(glmnet.obj$lambda==optimal.lambda)
# take beta for optimal lambda
optimal.beta - glmnet.obj$beta[,lambda.index]
# find non zero beta coef
nonzero.coef - abs(optimal.beta)0
selectedBeta - optimal.beta[nonzero.coef]
# take only covariates for which beta is not zero
selectedVar - predict_matrix[,nonzero.coef]
# create a dataframe for trainSet with time, status and selected
variables in binary representation for evaluation in pec
reformat_dataSet - as.data.frame(cbind(surv_obj,selectedVar))
# glmnet.cox only with meaningful features selected by stepwise
bidirectional AIC feature selection
glmnet.cox.meaningful - step(coxph(Surv(time,status) ~
.,data=reformat_dataSet),direction=both)
##--
-
##MODEL PERFORMANCE
##--
-
##
## Calculate the Brier score - pec does its own bootstrap so this
function runs on i=51 (i.e., whole trainset)
## Brier score calculation to cox-glmnet
peperr.coxglmnet - peperr(response=surv_obj, x=selectedVarCox,
fit.fun=fit.coxph, load.all=TRUE,
indices=resample.indices(n=nrow(surv_obj),
method=boot, sample.n=50))
# Get error predictions from peperr
prederr.coxglmnet - perr(peperr.coxglmnet)
# Integrated prediction error Brier score calculation
ipec.coxglmnet-ipec(prederr.coxglmnet,
eval.times=peperr.coxglmnet$attribute, response=surv_obj)
## C-Index calculation 50 iter bootstrapping
for (i in 1:50){
print (paste(Iteration:,i))
train - sample(1:nrow(dataset), nrow(dataset), replace =
TRUE) ##
random sampling with replacement
# create a dataframe for trainSet with time, status and
selected
variables in binary representation for evaluation in pec
reformat_trainSet - reformat_dataSet [train,]
# glmnet.cox
Re: [R] Interperting results of glmnet and coxph plot, Brier score and Harrel's C-Index - am I doing something wrong ???
This entire procedure is not valid. You cannot use a penalized method for
selecting variables then use an unpenalized procedure on those selected.
Frank
David Winsemius wrote
On Sep 28, 2013, at 2:39 AM, E Joffe wrote:
Hi all,
I am using COX LASSO (glmnet / coxnet) regression to analyze a
dataset of
394 obs. / 268 vars.
I use the following procedure:
1. Construct a coxnet on the entire dataset (by cv.glmnet)
2. Pick the significant features by selecting the non-zero coefficient
under the best lambda selected by the model
3. Build a coxph model with bi-directional stepwise feature selection
limited to the coxnet selected features.
I was a bit puzzled by the third step. Once you had a reduced model
from glmnet, what was the statistical basis for further elimination of
variables?
(Quite apart from the statistical issues, I was rather surprised that
this procedure even produced results since the 'step' function is not
described in the 'stats' package as applying to 'coxph' model objects.)
To validate the model I use both Brier score (library=peperr) and
Harrel's [Harrell]
C-Index (library=Hmisc) with a bootstrap of 50 iterations.
MY QUESTION : I am getting fairly good C-Index (0.76) and Brier
(0.07)
values for the models however per the coxnet the %Dev explained by
the model
is at best 0.27 and when I plot the survfit of the coxph the plotted
confidence interval is very large.
How many events did you have? (The width of CI's is most importantly
dependent on event counts and not particularly improved by a high case
count. The power considerations are very similar to those of a
binomial test.)
What am I missing here ?
Perhaps sufficient events? (You also seem to be missing a description
of the study goals.)
--
David.
%DEV=27%
Brier score - 0.07 ($ipec.coxglmnet - [1] 7.24)
C-Index - 0.76 ($cIndex - [1] 0.763)
DATA: [Private Health Information - can't publish] 394 obs./268 vars.
CODE (need to define a dataset with 'time' and 'status' variables):
library(survival)
library (glmnet)
library (c060)
library (peperr)
library (Hmisc)
#creat Y (survival matrix) for glmnet
surv_obj - Surv(dataset$time,dataset$status)
## tranform categorical variables into binary variables with
dummy for
dataset
predict_matrix - model.matrix(~ ., data=dataset,
contrasts.arg = lapply
(dataset[,sapply(dataset, is.factor)], contrasts))
## remove the statu/time variables from the predictor matrix (x)
for
glmnet
predict_matrix - subset (predict_matrix, select=c(-time,-status))
## create a glmnet cox object using lasso regularization and cross
validation
glmnet.cv - cv.glmnet (predict_matrix, surv_obj, family=cox)
## get the glmnet model on the full dataset
glmnet.obj - glmnet.cv$glmnet.fit
# find lambda index for the models with least partial likelihood
deviance (by cv.glmnet)
optimal.lambda - glmnet.cv$lambda.min# For a more parsimoneous
model use lambda.1se
lambda.index - which(glmnet.obj$lambda==optimal.lambda)
# take beta for optimal lambda
optimal.beta - glmnet.obj$beta[,lambda.index]
# find non zero beta coef
nonzero.coef - abs(optimal.beta)0
selectedBeta - optimal.beta[nonzero.coef]
# take only covariates for which beta is not zero
selectedVar - predict_matrix[,nonzero.coef]
# create a dataframe for trainSet with time, status and selected
variables in binary representation for evaluation in pec
reformat_dataSet - as.data.frame(cbind(surv_obj,selectedVar))
# glmnet.cox only with meaningful features selected by stepwise
bidirectional AIC feature selection
glmnet.cox.meaningful - step(coxph(Surv(time,status) ~
.,data=reformat_dataSet),direction=both)
##--
-
##MODEL PERFORMANCE
##--
-
##
## Calculate the Brier score - pec does its own bootstrap so this
function runs on i=51 (i.e., whole trainset)
## Brier score calculation to cox-glmnet
peperr.coxglmnet - peperr(response=surv_obj, x=selectedVarCox,
fit.fun=fit.coxph, load.all=TRUE,
indices=resample.indices(n=nrow(surv_obj),
method=boot, sample.n=50))
# Get error predictions from peperr
prederr.coxglmnet - perr(peperr.coxglmnet)
# Integrated prediction error Brier score calculation
ipec.coxglmnet-ipec(prederr.coxglmnet,
eval.times=peperr.coxglmnet$attribute, response=surv_obj)
## C-Index calculation 50 iter bootstrapping
for (i in 1:50){
print (paste(Iteration:,i))
train - sample(1:nrow(dataset), nrow(dataset),
