0
$\begingroup$

I am currently using genomic expression levels, age, and smoking intensity levels to predict the number of days Lung Cancer Patients have to live. I have a small amount of data; 173 patients and 20,438 variables, including gene expression levels (which make up for 20,436). I have split up my data into test and training, utilizing an 80:20 ratio. There are no missing values in the data.

I am using knn() to train the model. Here is what the code looks like:

prediction <- knn(train = trainData, test = testData, cl = trainAnswers, k=1)

Nothing seems out of the ordinary until you notice that k=1. "Why is k=1?" you may ask. The reason k=1 is because when k=1, the model is the most accurate. This makes no sense to me. There are quite a few concerns:

  1. I am using knn() to predict a continuous variable. I should be using something along the lines of, cox maybe.
  2. The model is waaaaaaay too accurate. Here are a few examples of the test answer and the model's predictions. For the first patient, the number of days to death is 274. The model predicts 268. For the second patient, test: 1147, prediction: 1135. 3rd, test: 354, prediction: 370. 4th, test: 995, prediction 995. How is this possible? Out of the entire test data, the model was only off by and average of 9.0625 days! The median difference was 7 days, and the mode was 6 days. Here is a graph of the results: Bar Graph.

So I guess my main question is what does knn() do, what does k represent, and how is the model so accurate when k=1? Here is my entire code (I am unable to attach the actual data):

# install.packages(c('caret', 'skimr', 'RANN', 'randomForest', 'fastAdaboost', 'gbm', 'xgboost', 'caretEnsemble', 'C50', 'earth'))
library(caret)

# Gather the data and store it in variables
LUAD <- read.csv('/Users/username/Documents/ClinicalData.csv')
geneData <- read.csv('/Users/username/Documents/GenomicExpressionLevelData.csv')
geneData <- data.frame(geneData)
row.names(geneData) = geneData$X
geneData <- geneData[2:514]
colNamesGeneData <- gsub(".","-",colnames(geneData),fixed = TRUE)
colnames(geneData) = colNamesGeneData

# Organize the data
# Important columns are 148 (smoking), 123 (OS Month, basically how many days old), and the gene data. And column 2 (barcode).
LUAD = data.frame(LUAD$patient, LUAD$TOBACCO_SMOKING_HISTORY_INDICATOR, LUAD$OS_MONTHS, LUAD$days_to_death)[complete.cases(data.frame(LUAD$patient, LUAD$TOBACCO_SMOKING_HISTORY_INDICATOR, LUAD$OS_MONTHS, LUAD$days_to_death)), ]
rownames(LUAD)=LUAD$LUAD.patient
LUAD <- LUAD[2:4]

# intersect(rownames(LUAD),colnames(geneData))
# ind=which(colnames(geneData)=="TCGA-778-7167-01A-11R-2066-07")
gene_expression=geneData[, rownames(LUAD)]

# Merge the two datasets to use the geneomic expression levels in your model
LUAD <- data.frame(LUAD,t(gene_expression))
LUAD.days_to_death <- LUAD[,3]
LUAD <- LUAD[,c(1:2,4:20438)]
LUAD <- data.frame(LUAD.days_to_death,LUAD)

set.seed(401)

# Number of Rows in the training data (createDataPartition(dataSet, percentForTraining, boolReturnAsList))
trainRowNum <- createDataPartition(LUAD$LUAD.days_to_death, p=0.8, list=FALSE)

# Training/Test Dataset
trainData <- LUAD[trainRowNum, ]
testData <- LUAD[-trainRowNum, ]

x = trainData[, c(2:20438)]
y = trainData$LUAD.days_to_death
v = testData[, c(2:20438)]
w = testData$LUAD.days_to_death

# Imputing missing values into the data
preProcess_missingdata_model <- preProcess(trainData, method='knnImpute')
library(RANN)
if (anyNA(trainData)) {
    trainData <- predict(preProcess_missingdata_model, newdata = trainData)
}
anyNA(trainData)

# Normalizing the data
preProcess_range_model <- preProcess(trainData, method='range')
trainData <- predict(preProcess_range_model, newdata = trainData)
trainData$LUAD.days_to_death <- y
apply(trainData[,1:20438], 2, FUN=function(x){c('min'=min(x), 'max'=max(x))})

preProcess_range_model_Test <- preProcess(testData, method='range')
testData <- predict(preProcess_range_model_Test, newdata = testData)
testData$LUAD.days_to_death <- w
apply(testData[,1:20438], 2, FUN=function(v){c('min'=min(v), 'max'=max(v))})

# To uncomment, select the text and press 'command' + 'shift' + 'c'
# set.seed(401)
# options(warn=-1)
# subsets <- c(1:10)
# ctrl <- rfeControl(functions = rfFuncs,
#                    method = "repeatedcv",
#                    repeats = 5,
#                    verbose = TRUE)
# lmProfile <- rfe(x=trainData[1:20437], y=trainAnswers,
#                  sizes = subsets,
#                  rfeControl = ctrl)
# lmProfile

trainAnswers <- trainData[,1]
testAnswers <- testData[,1]

library(class)
prediction <- knn(train = trainData, test = testData, cl = trainAnswers, k=1)

#install.packages("plotly")
library(plotly)
Test_Question_Number <- c(1:32)
prediction2 <- data.frame(prediction[1:32])
prediction2 <- as.numeric(as.vector(prediction2[c(1:32),]))
data <- data.frame(Test_Question_Number, prediction2, testAnswers)
names(data) <- c("Test Question Number","Prediction","Answer")

p <- plot_ly(data, x = ~Test_Question_Number, y = ~prediction2, type = 'bar', name = 'Prediction') %>%
    add_trace(y = ~testAnswers, name = 'Answer') %>%
    layout(yaxis = list(title = 'Days to Death'), barmode = 'group')
p
merge <- data.frame(prediction2,testAnswers)

difference <- abs((merge[,1])-(merge[,2]))
difference <- sort(difference)
meanDifference <- mean(difference)
medianDifference <- median(difference)
modeDifference <- names(table(difference))[table(difference)==max(table(difference))]
cat("Mean difference:", meanDifference, "\n")
cat("Median difference:", medianDifference, "\n")
cat("Mode difference:", modeDifference,"\n")

Lastly, for clarification purposes, ClinicalData.csv is the age, days to death, and smoking intensity data. The other .csv is the genomic expression data. The data above line 29 doesn't really matter, so you can just skip to the part of the code where it says "set.seed(401)".

Some samples of the data (Days to death, A number representative of age, and the rest are genes and their expression levels. I only shared a few. There are over 20,000):

days_to_death    OS_MONTHS
121              3.98

NACC1   2001.5708   2363.8063   1419.879
NACC2   58.2948     61.8157     43.4386
NADK    706.868     1053.4424   732.1562
NADSYN1 1628.7634   912.1034    638.6471
NAE1    832.8825    793.3014    689.7123
NAF1    140.3264    165.4858    186.355
NAGA    1523.3441   1524.4619   1858.9074
NAGK    983.6809    899.869     1168.2003
NAGLU   621.3457    510.9453    1172.511
NAGPA   346.9762    257.5654    275.5533
NAGS    460.7732    107.2116    321.9763
NAIF1   217.1219    202.5108    132.3054
NAIP    101.2305    87.8942     77.261
NALCN   13.9628     36.7031     48.0809
NAMPT   3245.6584   1257.8849   5465.6387

Edit: I just realized OS_Months is actually the days to death. It means Overall Survival Months, basically how many months the patient has to live. I am reformatting my code accordingly. I will re-edit this when I figure out whether it is the problem or not.

$\endgroup$
  • $\begingroup$ Do you have multiple clinical observations per patient? (eg. for the 1st patient, two observations that are 6 days apart?) $\endgroup$ – Matthew Gunn Jul 3 '18 at 18:28
  • $\begingroup$ Nope, it's 1 observation per patient. $\endgroup$ – Octavius Jul 3 '18 at 18:45
  • $\begingroup$ Possible causes of fantastical, impossible accuracy include (1) training on your test data or (2) using predictors that are essentially what you're trying to predict (eg. regress left shoe size on right shoe size). I'm concerned (like you are) that something like that is going on? $\endgroup$ – Matthew Gunn Jul 18 '18 at 15:30
1
$\begingroup$

To answer some of your questions

I guess my main question is what does knn() do, what does k represent, and how is the model so accurate when k=1?

knn() finds the k records in your dataset (the k-nearest neighbors) that are closest to the record it is currently trying to classify. What we mean by closest is that the distance between the records calculated using your auxiliary variables and some distance measure (knn probably defaults to Euclidian distance but I am not sure on that). It then takes the average of the k-nearest neighbors and that is what it uses to classify the record of interest.

K represents the number of records to include in your neighborhood. If you choose k=1 you look only at the record that is closest.

Typically bias will increase as the choice of K increases but variance will decrease.

As to why your model is so accurate, without looking at your data I would not be able to tell you. It could be that you simply have very predictive auxiliary information

$\endgroup$
  • $\begingroup$ I love the way you explain knn(). I finally understand it now. Maybe not Euclidian distance, but everything else makes sense. Thank you! $\endgroup$ – Octavius Jul 5 '18 at 13:16
  • $\begingroup$ By the way, I have posted a sample of my data at the end of my post; it's 1 patient. I was not able to post all of the patient's data because there would be 20,000 genes, so I posted like 5-7 genes. Anyways, it's basically like that for like, 500 patients. $\endgroup$ – Octavius Jul 5 '18 at 13:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.