Stochastic approaches play a vital role in weather, climate, and, more in general, geophysics systems, addressing processes and scales beyond the resolution of deterministic models. Similar to equilibrium/non-equilibrium thermodynamics, intricate fast and local dynamics may not always be the primary focus. Practical applications often prioritize observables capturing phenomena at dominant temporal and spatial scales. Developing models for these "large-scale" observables, resulting from averaging fast and local contributions, can be simplified into Low Order Models (LOMs) with reduced degrees of freedom described by ordinary differential equations. Unresolved degrees of freedom are introduced as stochastic components, exhibiting either Markovian or non-Markovian characteristics. The challenge lies in deriving dependable stochastic differential equations representing the statistics of real large-scale, slow features in the climate/ocean system. While paralleling material physics, it is crucial to recognize that direct transfer of tools and outcomes is hindered by the non-Hamiltonian nature of climate/geophysical LOMs and the impracticality of a Markovian treatment of noise due to wide-ranging time scales. A critical examination of the conventional statistical mechanics approach, customized for such LOMs, becomes essential. To this end, we propose utilizing an approach based on the operator cumulant method, which has been recently revisited and generalized, along with the linear response method in a non-Hamiltonian setting. Formal results are then derived, and applications to some typical classes of examples are presented to clarify this approach.