# How to properly smooth a 2d map?

Lately I've been trying to model the probability of success for an event given it's location. Basically, the data I have consists of locations in x and y coordinates, and a corresponding value of 1 or 0 (indicating success or failure). When I bin the events into boxes and calculate their average success rate for that bin, I get:

My goal is to smooth this out, as there is obviously some noise in here. I've tried multiple methods (nonlinear regression, logit model, mgcv package) of which none seemed to give a proper smooth. For example, a full tensor product smooth from the binomial family (using the mgcv package), yielded this result:

This is kind of in the right direction, but it does seem to underfit the very high probabilities around the 'white' area. Upping the degrees of freedom (currently 15) in the model is limited due to computatoinal complexity.

I've been trying to get a more accurate model of the probabilities for some weeks without success. Any tips/ideas on how to improve this? (I'm using R at the moment.)

• is this a point process? i.e., the number of points and their location are random? A realization of the process would be a set of points $\{(x_1, y_1), \ldots,(x_m, y_m)\}$. I'm not familiar with point processes, but it looks like your problem is that you want less smoothing, without increasing the degrees of freedom even more. Have you already tried local smoothing? See ?adaptive.smooth from mgcv. Other than this, for an answer I guess we would need more domain knowledge. What physical/social system are you modeling? Are events independent? Etc. Commented Nov 9, 2016 at 23:26
• Thanks for the tip, I'll take a look into local smoothing. What I'm modeling is the probability of a shot being a goal in football. The data is just shot locations and a 1 or 0, indicating if it was a goal or not. Shots can be assumed to be independent. Since players only shoot close to goal, the sample size for shots far away or very very close (which is difficult) is rather small. If you have any additional comments I'd love to hear about them. Commented Nov 10, 2016 at 13:58
• I'd love to help, but I know very little about football :) I'll try to see if I can help you get an answer from someone else. If I understand correctly, you have a finite grid of shot locations (even if the average soccer player is very small with respect to the size of the soccer field, I don't think it makes sense to distinguish between shot locations 1 cm apart...). Denoting with $(x_i,y_i)$ the center of grid cell $i$, you define a set of Bernoulli variables $\{X(x_1,y_1),\ldots,X(x_m,y_m)\}$, each one associated to the probability of scoring from cell $i$... Commented Nov 11, 2016 at 10:10
• ...you have an i.i.d. sample for each cell and you want to estimate the $m$ scoring probabilities. So you don't need to predict from where the player will shoot, but only whether he will score, given the location from which he's shooting. Right? Commented Nov 11, 2016 at 10:14
• I think you should look at beta regression with spatial predictors - have a look at this question. Commented Nov 11, 2016 at 10:27

I do wonder if there is a big problem here at all? If you look at the very white pixels in your upper figure, they are surrounded by points that are less intensely-coloured (?) white. What you have here are data and the fitted spline is balancing the highly successful point with the less successful ones in the same spatial location.

A basis dimensionality of 15 sounds far too few for a 2-d thin plate or tensor product spline. In this are you referring to the k in s(x, y, k = 15) or te(x, y, k = 15), or did you set k to something else? If the former, then you actually have 15*15 as the basis dimension, which would seem more reasonable as a starting point: the smoothness selection will reduce this via the wiggliness penalty.

As you have both known 0 and 1's, I would start, if using mgcv, using the actual coordinates of the shots (not aggregated to pixels) with:

m <- gam(scored ~ te(x, y), data = myDF, family = binomial, method = "REML")


and then check the resulting fit for adequacy, especially via gam.check(). This will tell you whether the basis dimension for the te() term was high enough or not (look for K close to 1 or a high p value in the printed output). You can also assess the model diagnostic plots but some of these may be less than useful as the response is binary.

If gam.check() suggests too low a basis dimension (the default is $2^5$) and you should increase k.

It might be that you need an adaptive smoother: you can do te(x, y, bs = "ad") but this will impact the computational efficiency a lot as this is effectively going to fit several different GAMs to section of the data.

If you are on Linux you could consider setting the nthreads option via gam.control to a number equal to the available CPU cores (or thereabouts) in the gam() call: control = gam.control(nthreads = 4). This will utilise multiple threads in some parts of the computation and can for some models greatly decrease the compute time.

You could also look at bam() as a drop-in replacement for gam(), which is designed to work on big data sets.