# Predicting a dichotomous variable

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 predict

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)


Which gives:

$fit 1 -2.774538$se.fit
[1] 1.492594

$residual.scale [1] 1  So... what do I do with this? ?predict.glm says $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!

-
I would first include only important variables and see if there is an interaction. See that the residuals of models are OK. Use AIC (or other tests) to see if you've improved your model by incorporating/removing the interaction(s). Once you have a valid model you can go on predicting. Make sure you use type = "response" in your predict to transform the predicted values from logit to odds. – Roman Luštrik Oct 17 '12 at 8:02
@Roman Luštrik: ahhhhh type=response! I did not know about that one! So that number I get is the logit, now it's starting to make much more sense – nico Oct 17 '12 at 9:16
You can transform predict values either by exp(x)/(1+exp(x)) (x is fitted values), plogis(model\$fit) or using the aforementioned type = "response". – Roman Luštrik Oct 17 '12 at 16:07
@FrankHarrell, can you shortly expand the reasoning behind not deleting "insignificant" variables? – Roman Luštrik Oct 18 '12 at 7:02
There are at least three safe approaches as described in my text Regression Modeling Strategies or my course notes which are at biostat.mc.vanderbilt.edu/rms. (1) Fit a very large model using penalized maximum likelihood estimation. (2) Use data reduction blinded to Y to reduce the number of predictors. (3) Use subject matter expertise to reduce the number of variables to what the effective sample size will support. – Frank Harrell Oct 21 '12 at 22:24