Best method to predict binary outcome with multiple records per subject I am interested in building a model to predict the binary outcome, retention (1 - retained; 0 - not retained) with various potential predictor variables (either continuous or categorical). 
With that being said, I have a dataset containing multiple records (magazine subscriptions) for some subjects. For example, I have 4 records at the magazine subscription level for Joe Smith (3 of which he is retained, while the other record he is not retained) and 7 magazine subscription records for John Doe (4 of which he is retained, while the other 3 records he is not retained). 
My initial thoughts were to use logistic regression using a single randomly selected magazine subscription for each subject. For example, randomly selecting 1 of the 4 magazine subscriptions for Joe Smith and 1 of the 7 magazine subscriptions for John Doe. Obviously, I would be losing a great deal of my data, which makes me think that there has to be a better way of modeling this data.
What method would you be best to predict retention with data such as this?
Thank you!
 A: I'm not an expert on this but to me it looks like you should use a logistic mixed effects regression model (i.e., with a logit link). This lets you declare a random intercept (and slope if required) for every subject and thus account for the dependence of the observations. This way you can use all the records and preserve all the data and variance. In R, lmer in the lme4 package can do this by choosing the "binomial" family. The model would look roughly like this (if using only a random intercept):
library(lme4)
model <- lmer(retention <- predictor1 + predictor2 + (1|subjectid), family=binomial)

I hope this is what you were looking for.
Edit: Because I'm currently working on a similar problem, I'll have to update this with some current insights. It is probably not a good idea to just modell a random intercept by (+ 1|subject) but also to assume that your one or all of your predictor variables may have a different effect for every subject by putting (1+predictor1+predictor2|subjectid), thus additionally modelling a random slope.
A: An almost optimal approach can be to fit ordinary binary regression model, possibly with updated (time-dependent) covariates, to all the records, and to just not trust the standard errors that result.  You can use the cluster sandwich covariance matrix estimator or cluster bootstrap to get "honest" standard errors.  I say "almost" because a mixed effects or other full likelihood approach that takes the particular form of intra-subject correlation into account can result in more efficient estimates of $\beta$ coefficients.  But this "working independence model" approach is often not far from optimal, as long as you compute the variance-covariance matrix so as to take into account redundancies across records within subject.  In the R rms package the pertinent functions are lrm, robcov, bootcov.
A: A large part of the problem seems to be that the data structure is not amenable to analysis. Is there any reason you can't collapse the data into one record per person-year (e.g you're using more than just the variables from the subrecords)? The result may have a heck of a lot of indicator variables, but it would make the analysis a lot easier.
Due to the different requirements of the data at various stages at my company, I've had to collapse data in that manner countless times.
