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.
