I have a series of descriptors, some continuous, some discrete and an output variable which is dichotomous.
I have several parameters, but for the sake of simplicity let's say my data look like:
Sex | Age | Genotype | Dose 1 | Dose 2 | Outcome ------|-------|------------|---------|-----------|------------ M | 32 | AABB | 150 | 30 | YES F | 65 | aaBb | 110 | 30 | YES M | 42 | AaBb | 200 | 50 | NO ...
I would like to make a predictive model to determine the optimal combination of
Dose 1 and
Dose 2 to have a desired outcome.
So my question is, if I have a new male subject of given genotype and age, what is the best combination of doses that will give a positive outcome with the highest probability? Or, to see things the other way around, given the other parameters, what are the odds of having a positive outcome with a given set of doses?
I thought I could use R to generate a linear model with
glm, and then use
predict to predict the outcome. However, I never really dealt with this type of problems before and I am a bit confused on how to use these functions and interpret their results. Also, I am not sure if this is the correct way to deal with the problem.
For instance, let's generate some random data:
set.seed(12345) num.obs <- 50 sex <- sample(c("M", "F"), num.obs, replace=T) age <- sample(20:80, num.obs, replace=T) genotype <- sample(LETTERS[1:8], num.obs, replace=T) dose.1 <- sample(100:200, num.obs, replace=T) dose.2 <- sample(30:70, num.obs, replace=T) outcome <- sample(0:1, num.obs, replace=T) data <- data.frame(sex=sex, age=age, genotype=genotype, dose.1=dose.1, dose.2=dose.2, outcome=outcome)
Which gives 50 observation such as
> head(data) sex age genotype dose.1 dose.2 outcome 1 F 78 C 183 54 0 2 F 70 E 156 66 1 3 F 39 H 180 35 0 4 F 32 E 135 51 0 5 M 64 E 121 57 1 6 M 50 H 179 61 1
Now, I generate a model with
model <- glm(outcome ~ sex + age + genotype + dose.1 + dose.2, data=data, family="binomial")
First question: without any a priori knowledge of the interactions between the descriptors, how do I choose the correct formula? Should I try various interactions and see which models gives the best fit e.g. looking at residual deviance or AIC? Are there functions to do this for me or should I try all of the combinations manually?
OK, let's say I found the model is good, now I use
new.data <- list(sex=factor("M", levels=c("M", "F")), age=35, genotype=factor("C", levels=LETTERS[1:8]), dose.1=150, dose.2=30) outcome <- predict(model, new.data, se=T)
$fit 1 -2.774538 $se.fit  1.492594 $residual.scale  1
So... what do I do with this?
$fit is the prediction but obviously that is not a yes/no type of prediction... what I would ideally need is something on the lines of "89% YES / 11% NO".
How do I interpret the result of
predict and how would I go about having the type of result I want?
Finally, are there functions to explore the parameter space so that I get a graph with the outcome in the dose1 vs dose2 space?
EDIT: just to clarify: I do not have a specific reason to use a generalized linear model, it is just something that came to my mind as a possible solution. If anyone can think of other ways to solve this type of problem I would gladly welcome their suggestion!