We study maximal operators associated to singular averages along finite subsets $Sigma$ of the Grassmannian $mathrm{Gr}(d,n)$ of $d$-dimensional subspaces of $mathbb R^n$. The well studied $d=1$ case corresponds to the the directional maximal function with respect to arbitrary finite subsets of $mathrm{Gr}(1,n)=mathbb S^{n-1}$. We provide a systematic study of all cases $1leq d<n$ and prove essentially sharp $L^2(mathbb R^n)$ bounds for the maximal subspace averaging operator in terms of the cardinality of $Sigma$, with no assumption on the structure of $Sigma$. In the codimension $1$ case, that is $n=d+1$, we prove the precise critical weak $(2,2)$-bound. Drawing on the analogy between maximal subspace averages and $(d,n)$-Nikodym maximal averages, we also formulate the appropriate maximal Nikodym conjecture for general $1<d<n$ by providing examples that determine the critical $L^p$-space for the $(d,n)$-Nikodym problem. Unlike the $d=1$ case, the maximal Kakeya and Nikodym problems are shown not to be equivalent when $d>1$. In this context, we prove the best possible $L^2(mathbb R^n)$-bound for the $(d,n)$-Nikodym maximal function for all combinations of dimension and codimension. Our estimates rely on Fourier analytic almost orthogonality principles, combined with polynomial partitioning, but we also use spatial analysis based on the precise calculation of intersections of $d$-dimensional plates in $mathbb R^n$.