# Smoothing a surface with missing and non-uniformly weighted data

I have an uneven data set as shown in the image below containing x (2 thru 7.5) and y (0 to 200) predictors. The circles represent values from 0 (light) to 1 (dark). I also have a value n for each point, representing the number of observations. This should likely be used to weight the samples. My csv only has a record for each circle and is not therefore a matrix. I'd like to use an R method to create a smooth "surface" where I can input an x and y and return a z from 0 to 1.

I have looked at the following question: https://stackoverflow.com/questions/20848740/smoothing-surface-plot-from-matrix My sense is this won't work given the absence of some data.

Your sense is not quite right; there is nothing stopping you using the data you show above with the GAM approach I discuss in the SO answer you link to.

Here you'll certainly want to use te(x, y) as a term in the model formula (not s(x, y)) because your x and y are scaled quite differently and we probably don't want to assume the same degree of wiggliness in both dimensions (variables).

You can weight the observations in the manner you suggest via the weights argument. Say you have the number of observations per x, y pair in variable n, then you can add weights = n to account for the number of observations per data point.

For the continuous response bounded at 0 and 1, you could consider this as a model for a response variable that is the proportion of dark as a function of (x, y) covariates. In that case a beta response model may be appropriate. If so, or to give it a try, you can add family = betar to the gam() call: see ?mgcv::betar for more details. (You should also use method = "REML" or method = "ML" in the gam() call - you may be required to with this extended GAMs.)

You call therefore is going to look something like

gam(z ~ te(x, y), data = myDF, family = betar(link = "logit"), method = "REML", weights = n)


### Some comments on the uneven distribution of data in x and y

What mgcv will do here is create a marginal basis for x and for y and use a tensor product of these bases to give the final 2d basis used in fitting the model. In this example, this may be a little wasteful of basis functions as for some combinations and x and y there is no support from the data, yet the tensor product of the marginal bases will cover the empty area of the data. There, some basis functions will evaluate to zero for all observations and hence are redundant and it would be useful to remove them before model fitting so as to reduce the number of parameters that need to be estimated.

Simon Wood, author of mgcv, has recently published a paper (Wood, 2016) discussing this issue, the results of which you can I believe take advantage of in the B spline basis available in newer versions of the package. The key result is, for simulated data, that the tensor product of marginal B spline basis uses far fewer coefficients than the full rank basis whilst achieving essentially the same fit as the full rank version:

The image (from Wood, 2016) shows the default full rank tensor product spline on the right and the reduced rank version on the left for the simulated data. The fit on the right required estimation of 625 coefficients, whilst the one on the left only 328.

Wood, S. N. P-splines with derivative based penalties and tensor product smoothing of unevenly distributed data. Stat. Comput. 1–5 (2016). doi: 10.1007/s11222-016-9666-x

The paper is available under an CC-BY open access licence so you can check it out if this interests you.