Methods of examining the fit of multi-dimensional point process models using residual analysis are proposed. One method involves rescaled residuals, obtained by transforming points along one coordinate to form a homogeneous Poisson process inside a random irregular boundary. Both vertical and horizontal forms of this rescaling are discussed. We also present a different method of residual analysis, involving thinning the point process according to the conditional intensity to form a homogeneous Poisson process on the original, untransformed space. These methods for assessing goodness-of-fit are applied to point process models for the space-time-magnitude distribution of earthquake occurrences, using in particular the multi-dimensional versio of Ogata's epidemic-type aftershock sequence (ETAS) model and a 30-year catalog of 580 earthquakes occurring in Bear Valley, California, as an example. The thinned residuals suggest that the fit of the model may be significantly improved by using an anisotropic spatial distance function in the estimation of the spatially varying background rate. Using rescaled residuals, it is shown that the temporal-magnitude distribution of aftershock activity is not separable, and that in particular, in contrast to the ETAS model, the triggering density of earthquakes appears to depend on the magnitude of the secondary events. The residual analysis highlights that the fit of the space-time ETAS model may be substantially improved by allowing the parameters governing the triggering density to vary for earthquakes of different magnitudes. Such modifications are important since the ETAS model is widely used seismology for hazard analysis.