A regression model whose response variable is the day of year that an annual event (usually) occurs In this particular case I'm referring to the day on which a lake freezes. This "ice-on" date only occurs once a year, but sometimes it doesn't occur at all (if the winter is warm). So on one year the lake may freeze on day 20 (january 20th), and another year it might not freeze at all.
The goal is to figure out drivers of ice-on date.
Predictors would be things like fall/ winter air temperature each year. Year could be a predictor for the long-term linear trend.
1) Is the integer "day of year" a reasonable response variable (if not, what is?)?
2) How should one handle the years when the lake never froze?
Edit:
I don't know what the etiquette is here, but I figured I'd post the outcome of the suggestions I received. Here's the paper, open access. I got good feedback on the approach used, thanks @pedrofigueira and @cboettig . Of course, errors are my own.
 A: I think one can consider "day of the year" as a response variable to a multivariate regression. In order to handle years when the lake never froze I would simply consider that the day of freezing is larger than an observable lower limit which corresponds, for instance, to the day when ice content starts to melt (or melts completely, if you want to be very conservative). Theoretically it should freeze after that, or can freeze after that, but we do not know. This way you could use the data you collected on the different parameters to understand how the freezing day depends on them, if it was allowed to be later than the latest observable date. You can then use a Tobit model to handle simultaneously freezing days (corresponding to "normal" datapoints) and lower limits (corresponding to limits and thus a censored regression). 
In order to correctly include the measured lower limits in the analysis, you can use a censored regression model in which the dependent variable has a cut-off at the value of the lower limit. The above-mentioned Tobit model is appropriate for this case; it assumes the existence of an unobservable (latent) dependent variable $y_i^*$ which in our case corresponds to the freezing date if the winter extended indefinitely. The observable dependent variable $y_i$ (i.e. the measured lower limit on freezing date) is then taken to be equal to the latent variable in the absence of a lower limit $L_i$, and equal to the lower limit otherwise
\begin{eqnarray}
   y_i = \left\{
     \begin{array}{ll}
       y_i^* & \quad \mathrm{if} \quad \bar{\exists}\,L_i \:\: (\textrm{i.e.} \, y_i^* < L_i) \\
       L_i   & \quad \mathrm{if} \quad y_i^*\geq L_i
     \end{array}
   \right.
\end{eqnarray} 
The application of the Tobit model to handle observation-by-observation censoring, results in a log-likelihood function of the form
\begin{eqnarray}
\mathcal{L} = \sum_{i \,\in\, y_i^* < L_i} ln \left[ \phi\left(\frac{y_i-X_{ij}\beta_j}{\sigma}\right)/\sigma \right] \,+\, 
              \sum_{i \,\in\, y_i^*\geq L_i}ln \left[ \Phi\left(\frac{L_i-X_{ij}\beta_j}{\sigma}\right) \right] \, \,
\end{eqnarray}
where $\phi(.)$ and $\Phi(.)$ denote the probability and cumulative density functions, respectively, of the the standard normal distribution. The index $i$ runs on the observations and $j$ on the independent variables. The solution to the linear regression is the set of parameters $\beta_j$(including intercept) that maximizes the log-likelihood function. 
A: Day of year is one sensible predictor variable, and for that I think it is sensible to treat it as @pedrofigueira suggests.  
For other predictor variables you may need to be careful about how you represent time.  For instance, imagine you have air temperatures by day -- how would you model air temperature as a predictor of ice-on day?  I don't think comparing same day-of-year samples is sufficient. 
In any such analysis, I think it helps to write down what you think a plausible generating model (or models) of the data might be, (where some physics might be available as a guide).  For instance, a reasonable model might be to integrate the number of days below freezing, and when that integral passes a threshold (e.g. related to the thermal mass of the lake), ice-on occurs.  From such a model you can then ask what is a reasonable approximation and what isn't. 
For instance, day-of-year as predictor matters to that model only in so much as day of year is a good predictor of temperature.  Thus knowing only the day of the year, one would just have an average day-of-year corresponding to the ice-on threshold, with perhaps some normal distribution about it resulting from interannual temperature variations, and therefore looking for a trend in day-of-year is completely justified.
But if you know other variables like air-temp by day, you probably face dealing with somewhat more complicated model more directly. If you are just using the annual values (minimums? means?) than variable as a predictor of ice-on day also seems reasonable (by the same argument as above).   
A: For this problem you need two response variables.  One Boolean response that indicates whether the lake froze or not, and one integer response giving the day of the year, conditional on the indicator being true.  In years when the lake froze, both the Boolean and integer are observed.  In years when the lake didn't freeze, the Boolean is observed and the integer is not. You can use a logistic regression for the Boolean. The regression for day of the year could be an ordinary linear regression.  
The circular nature of the day of the year shouldn't be a problem as long as you number the possible freeze-over days consecutively within a given time period.  If you are wondering where to start the numbering, I would suggest the day when the predictors were measured.  If you want the model to represent causal effects, it must be the case that all predictors were measured before any possible freeze-over.
To handle the integer and bounded nature of the day of the year, could use a discretization model.  That is, there is a real latent value which generates an observation in the following way:  if the value is within the bounds then the observation equals the latent value rounded to the nearest integer, otherwise the value is truncated to the bounds.  The latent value itself can then be modelled as a linear function of the predictors plus noise.
A: What you have is time-to-event data, which is also termed survival analysis.  That is not really my area, so I am not giving a detailed answer here.  Googling for "time-to event data" or "survival analysis" will give you a lot of hits!
One good starting point could be the chapter (13) about survival analysis in Venables/Ripley: MASS, or the classic
"The Statistical Analysis of Failure Time Data, Second Edition"  by
John D. Kalbfleisch, Ross L. Prentice(auth.)
EDIT, EXTENDED ANSWER
As an alternative to survival analysis, you could approximate that by ordinal logistic regression.  By example, in your example case of first freezing date, define some dates for which you give the "have been freezing at or before" state, 0 (no freezing), 1 (freezing).
That nicely accomodates the years without freezing, you simply have an all-zero response vector.  If your choosen dates are, say, 
1:08   15:08 1:09 15:09 1:10 15:10 1:11 15:11 1:12  15:12  1:01  15:01
and the actual date of first freezing was  17:11, then your observed vector will be
0       0    0    0     0    0     0    0      1     1     1      1

and, in general, all response vectors will have an initial block of zeros, followed by a block of ones.  Then, you can use this with ordinal logistic regression, obtaining an
estimated probability of freezing for each date.  Plotting that curve will give an approximation for a survival curve (survival, in this context, becomes "not having frozen  yet").
EDIT

One could also see your data as recurrent events, since river freezes (almost) every year.  Se my answer here:  Finding significant predictors of psychiatric readmissions
