We give a hyperpfaffian formulation for correlation functions in $β$-ensembles of $M \times M$ random matrices when $β= L^2$ is an even square integer. More specifically, to the $m$th correlation function $R_m : \R^m \rightarrow [0, \infty)$ we associate the $L$-vector valued function $ω_m : \R^m \rightarrow Λ^L \R^{L(M-m)}$ such that $R_m(\mathbf y)$ is given by the Vandermonde determinant in $y_1, \ldots, y_M$ times the hyperpfaffian of $ω_m.$ The partition function of the ensemble was previously shown to be the hyperpfaffian of a {\it Gram} $L$-form $ω$ in $Λ^L \R^{LM},$ and we demonstrate the relationship between $ω_m(\mathbf y)$ and $ω$, both having coefficients built from integrals of Wronskians of monic polynomials. Assuming the existence of families of polynomials sympathetic with the weight of the ensemble, we may construct $ω(\mathbf y)$ so it is very sparse (relative to the expected ${L(M-m) \choose L}$ coefficients of a general $L$-vector). These generalize skew-orthogonal polynomials arising in the well-understood $β= 4$ situation. Finally we explore the situation in the circular $β= L^2$ ensembles. Here the monomials give a prototype, and we give explicit formulas for (the circular versions of) $ω$ and $ω_m.$ We use our hyperpfaffian framework to produce exact formulas for the two point function when $β= 16$ for small values $M.$ Along the way we will record hyperpfaffian evaluations using known values of partition functions of $β$-ensembles.