Fitting Beta Distribution Parameters to Y conditional on X

I have a bi-variate data set where Y is in [0,1]. X is some measure of intensity and in this example is (0,~200) though there is no hard upper bound. X has a strong positive skew but I am not interested in the distribution of X right now. I don't believe what I want is a bivariate distribution.

I'm looking for a way to fit Beta Parameters conditional on X. An additional constraint would be that the mean(Y|X) should be non-decreasing as X increases.

I've tried a rolling window over X and fitting Beta to each subset of Y. Finally running a Linear Model of Shape1 ~ X, and Shape2 ~ X.

I would like to know if I'm on the right path, or if there is a better way. My preferred platform is R, but any help is appreciated.

Beta regression as proposed by Ferrari & Cribari-Neto (2004, Journal of Applied Statistics) can be used to model the parameters of a beta regression by linear predictors. Specifically, the beta distribution is then not characterized by two shape parameters $$p$$ and $$q$$ but instead by mean $$\mu$$ and precision parameter $$\phi$$ where $$\mu = p/(p + q)$$ and $$\phi = p + q$$. These parameters are then linked to linear preditors by suitable link functions, e.g., $$\begin{eqnarray*} \mathrm{logit}(\mu_i) & = & x_i^\top \beta\\ \log (\phi_i) & = & x_i^\top \gamma \end{eqnarray*}$$ This model is implemented in the R package betareg (Cribari-Neto & Zeileis, Journal of Statistical Software, doi:10.18637/jss.v034.i02). However, this does not implement any monotonicity restrictions on the regression coefficients.