Best way to model a regression with proportions as dependent variable? Consider a simple model for school graduation: there is school $i = 1,2,..., N$, each getting $y_i$ of $n_i$ students graduated. Assume $y_i\sim Binomial(n_i, \theta_i)$ where $\theta_i$ is the probability of graduation of each student at school $i$.
What is the best way to model the relationship between the individual probability of graduating $\theta_i$ and the observed characteristics $X_i$ for each school?
One possibility is to model the probability $\theta_i = X_i\beta$. Are there better ways?
 A: Proportion is the mean value of the variable where 1 represents success and 0 represents failure. So I suspect data has already been aggregated to give you proportions. For example you might have information about school graduation proportion given predictor variables such as average school grade, average aptitude test results of the school and so on. So the best way would be to use logistic regression on non-aggregated  data (predicting graduation given the pupil grades and aptitude tests). You could model the effect of higher level units (schools) by using hierarchical logistic regression. 
In any case, I don't see any justification to model probability itself. It really doesn't sound statistically (or philosophically). 
Answer based on a question edit:
If you only have school level data and not individual data then you don't need hierarchical regression. Also logistic regression is a binomial regression with logistic link function (in the framework of generalized linear models). So you can implement your problem in the following way:
glm(proportion ~ x1+x2, family=binomial, weights=w)

Proportion is simply $y_i/n_i$. Weights are $n_i$ and $x_{1},x_{2}$ are explanatory variables on a school level. 
To clarify even further I will provide a working example. Suppose we test an insecticide effectiveness dependent on dose and sex of an insect. There are 20 insects in every experimental batch:
dose <- rep(0:5, 2)

sex <- factor(rep(c("M", "F"), c(6, 6)))

dead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)

We calculate mortality dependent on dose and sex of the insect with interaction term:
glm(dead/20 ~ sex*ldose,family=binomial, weights=rep(20,12))

Example works flawlessly. But if we specify equivalent set of weights in the following way:
glm(dead/20 ~ sex*ldose, family=binomial, weights=rep(19, 12))

Package provides an error message: non-integer #successes in a binomial glm! But you will see that results in both cases are the same as it should be. This is because relative weights in both cases are equivalent. Program gives a prudent warning since weights are usually $n_i$ and proportion is calculated as $y_i/n_i$. So if $n_iy_i/n_i $ is not an integer program produces a warning. But if you know what you are doing you can simply ignore the warning.
I hope that solves your problem.
