Proportion Estimates - Shrink to the Mean Based on Sample Size Set-up:
One of the fundamental elements of direct marketing is choosing which lists of prospects to select and send an offer to. There are hundreds of lists on the market to choose from and each list will contain a specific number of records you can send an offer to. 
I have a table of historical performance of list. For each, the number of times a promotion was sent to someone on that list, along with the number of responders is given. So, the table might look like this - it is a made up example with only a couple of lists versus the hundred actually but the very small response rates (p) are realistic.
lists<-data.frame(listName=c('A','B','C','D','E'),numberPromoted=c(1000,5000,25000,35000,1654),numberRespond=c(3,45,98,350,5))
lists$p<-lists$numberRespond/lists$numberPromoted

I am looking to use linear programming to optimally allocate list selection in future campaigns. So, given the current number of records available in each list at the preesent time and the best estimate of the response rate of the list (p), choose how many records from a list to select, given constraints.
Problem: 
My question is how to best go about estimating p. Given these small response rates, I would like to shrink those without "adequate" information to a mean value. The naive way is to say that there has to be at least X number of records promoted before we trust the p estimated empirically. I am wondering if there is any other method that makes sense? 
For example, does it make sense to treat the list as a random effect and shrink the estimate? I THINK this code is shrinking the estimates of the response rate (p) to the mean rate for smaller cells.
install.packages("lme4")
library(lme4)
install.packages("arm")
library(arm)

#make individual 1/0 records
Lists<-data.frame(
listName<-c(replicate(1000,"List A"),replicate(5000,"List B"),replicate(25000,"List C"),replicate(35000,"List D"),replicate(1654,"List E")),
outcome<-c(replicate(997,0),replicate(3,1),replicate(4955,0),replicate(45,1),replicate(24902,0),replicate(98,1),replicate(34650,0),replicate(350,1),
    replicate(1649,0),replicate(5,1))
)


    #random intercept model
    mod<-lmer(outcome~1 + (1|listName),data=Lists,family=binomial)

    invlogit(coef(mod)$listName[1,1])
    invlogit(coef(mod)$listName[2,1])
    invlogit(coef(mod)$listName[3,1])
    invlogit(coef(mod)$listName[4,1])
    invlogit(coef(mod)$listName[5,1])

List    Modeled P   "TRUE" P
A   0.004241    0.003000
B   0.008602    0.009000
C   0.003981    0.003920
D   0.009924    0.010000
E   0.003961    0.003023





  #compare raw versus shrunk modeled estimates

raw<-sqldf("select listName, avg(outcome) as raw_p from Lists group by listName")
    meanRaw<-sqldf("select avg(outcome) as total_p from Lists")
    modeled<-c(
    invlogit(coef(mod)$listName[1,1]),
    invlogit(coef(mod)$listName[2,1]),
    invlogit(coef(mod)$listName[3,1]),
    invlogit(coef(mod)$listName[4,1]),
    invlogit(coef(mod)$listName[5,1]))

    plot(x=as.numeric(raw$listName),y=raw$raw_p, col="blue")
    points(x=raw$listName,y=modeled, col="red")
    abline(h=meanRaw)

 A: This looks about right as a strategy.  There is presumably a population of lists of which these are a sample.  If so, they're a natural random effect candidate.  And you are currently shrinking to the mean.
You might even know something about the lists on, e.g. their source or the demographic they target or some such.  Then you could shrink towards group means for groups or features actually relevant to your problem, rather than just the grand mean.  If so, put that information on the right hand side of the bar instead of the list indicator.
Also, since the counts are very small you also might want to look at using a Poisson GLM instead of the binomial.  For that, predict the numberRespond with family=poisson and use numberPromoted as offset.  Then recover the proportions with a simple division.  In any case, I would recommend avoiding any applications of linear regression here.
On a computational note: you don't have to make individual records.  Binomial GLMs in base R's glm and in lme4 will take the aggregated output and give you the same results.
Finally, it is possible that your subsequent optimisation routines can be made to make use of your uncertainty about the estimated proportions.  If there is much you will probably do better by pushing it through to the later analysis.  In your current formulation ranef or a posterior sample will give you this information.
A: Another approach would be to "manually" shrink the response rate to the sample mean using a quadratic penalty for deviating from the mean. If one uses quadratic loss function (with weighting each list by numberPromoted to take into account the amount of data in each list) then this becomes a modified version of weighted ridge regression with shrinking to the mean as opposed to shrinking to zero (aka "poor man's bayesian")
The penalty coefficient can be tuned using cross-validation.
