I have a dataframe (merged_COIN_plink) with columns of 181 RSID numbers such as (rs2807367_G) - (example just shown with two of them).

I am trying to do a coxph using the survival package on Rstudio. I have managed to do this for each RSID individually but with 181 i wondered if there was a way to sort of parallelise it or run it so it automatically goes through each RSID.

This is the code for the individual coxph:

coxph(Surv(merged_COIN_plink$SURVIVAL, merged_COIN_plink$EVENT) ~ rs2807367_G PS RS WBC, data= merged_COIN_plink)

I've looked through other posts, but they all do seem confusing. I wondered if it would be possible to use a wildcard such thing like RS* but I'm not sure if that is easy to do on Rstudio. I also thought that on unix you can make a list of the RS numbers and run the code through that to pick each RSID out but I don't know if that is possible on Rstudio.

I tried to make a list out of the column names starting with RS but this didn't seem to work properly:

rs_list <- merged_COIN_plink[grep("^rs",colnames(merged_COIN_plink)),] View(rs_list)

I also wasn't sure if a for loop would work but couldn't figure it out for column names not named all the same.

structure(list(ID = c("100002", "100003", "100004", "100005", 
"100006", "100007", "100008", "100010", "100011", "100012", "100013", 
"100014", "10004", "1002", "1003", "1004", "1005", "1006", "1007", 
"1008", "1010", "101001", "101002", "101003", "101004"), PHENOTYPE = c(-9L, 
-9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, 
-9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L, -9L), rs2807367_G = c(1L, 
0L, 2L, 2L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 
0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L), rs34963268_C = c(1L, 1L, 1L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 0L, 1L, 0L, 2L, 
0L, 0L, 1L, 0L, 0L, 0L), EVENT = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 
1L, 1L), SURVIVAL = c(349L, 384L, 283L, 671L, 674L, 285L, 224L, 
687L, 571L, 495L, 510L, 302L, 159L, 44L, 85L, 347L, 604L, 447L, 
1230L, 444L, 1260L, 758L, 392L, 379L, 188L), PS = c(2L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 2L, 1L, 1L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 1L), RS = c(1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 
0L), WBC = c(8.8, 8.1, 9.3, 8.9, 7.2, 6.7, 11.6, 10.7, 6.1, 12.9, 
10.1, 9.1, 6.8, 13.3, 13.5, 10.9, 8.7, 11.4, 9.8, 8.9, 8, 11.3, 
6, 5.6, 8.8)), row.names = c(NA, 25L), class = "data.frame")

CodePudding user response：

Iterating over a for loop is certainly valid, but with so many columns in your (real) data may take a long time. An answer that may be a a bit more "elegant" :)

# Get names of columns with RSID numbers (could also use `grep` or some other way
RSIDcols <- colnames(merged_COIN_plink[,3:4]) 

# Define formulas
formulas <- sapply(RSIDcols, 
                   function(x) as.formula(paste('Surv(SURVIVAL, EVENT) ~ ', paste(x, "  PS   RS   WBC"))))

# Run models
models <- lapply(formulas, function(x) {coxph(x, data = merged_COIN_plink)})

This outputs each result in a list (models) which can then be accessed by RSID name with the following:

names(models) <- RSIDcols # rename list elements to RSID name

# Access individual model results using `[[]]`
models[["rs34963268_C"]]

# > models[["rs34963268_C"]]
# Call:
#   coxph(formula = x, data = merged_COIN_plink)
# 
# coef exp(coef) se(coef)      z      p
# rs34963268_C -0.11518   0.89120  0.36030 -0.320 0.7492
# PS            0.72277   2.06012  0.43391  1.666 0.0958
# RS           -1.02832   0.35761  0.54250 -1.896 0.0580
# WBC           0.05065   1.05195  0.11965  0.423 0.6721
# 
# Likelihood ratio test=6.04  on 4 df, p=0.1962
# n= 25, number of events= 21 

CodePudding user response：

A solution using a for loop which should be sufficient for what you need, as you don't have too many columns to loop over and speed seems unlikely to be an issue. I also think it's quite easy to see what's going on.

This way stores each model in a list.

# identify the columns you need 
columns = colnames(dat)[grep("rs", colnames(dat))]

# create a list to store the results in
res = list()

# loop over the columns and perform the regression and save results in the list
for (i in seq_along(columns)) {
    res[[i]] = coxph(Surv(SURVIVAL, EVENT) ~ get(columns[i])    PS   RS   WBC, data = dat)
}

You can then make a table of whatever values you like as follows:

data.frame(nevent = sapply(res, function(x) x$nevent), columns)

  nevent      columns
1     21  rs2807367_G
2     21 rs34963268_C