Predicting a dichotomous variable

Question

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!