Spatial Sampling Designs - Sampling Spatial Units for Agricultural Surveys - page 171

using the clara function from the cluster package. Then, we use the gam
function of the mgcv package to estimate a spline response surface (see Wood 2006
for a review of the mgcv package). Before selecting the cube, we apply the
matrix package to easily build a block diagonal matrix.

> kk < - 10 # number of knots
> lambda_sq < - 3 # penalty parameter
> set.seed(200694)
> library(cluster)
> knots < - clara(X, k ¼ kk)$medoids
> library(mgcv)
> mod < - gam(yobs ~ s(xc, yc, k ¼ kk), data ¼ framepop,
+ knots ¼ list(x ¼ knots[,1], y ¼ knots[,2]))
> Z < - model.matrix(mod)[,-1]
> X < - cbind(rep(1, dim(framepop)[1]), framepop$xc, framepop$yc,
+ framepop$xc^2, framepop$yc^2)
> Px < - X%*%solve(t(X)%*%X)%*%t(X)
> I < - diag(rep(1, dim(framepop)[1]))
> C < - cbind(X, (I-Px)%*%Z)
> Q < - diag(rep(1, dim(Z)[2]))
> A1DA2 < - svd(solve(t(Z)%*%(I-Px)%*%Z+
+ lambda_sq*solve(Q))%*%t(Z)%*%(I-Px)%*%Z)
> A1DA2$d
[1] 9.970179e-01 9.969821e-01 9.968819e-01 9.756851e-01 9.752351e-01
9.016462e-01 8.917957e-01 1.720603e-12 4.529269e-13
> A1D < - A1DA2$u%*%diag(A1DA2$d)
> library(Matrix)
> B1 < - C%*%bdiag(diag(rep(1, dim(X)[2])), A1D)
> pik < - rep(n/N,N)
> B1 < - as.matrix(cbind(pik,as.matrix(B1)))
> bal3 < - samplecube(B1,pik,comment ¼ TRUE,method ¼ 1)

BEGINNING OF THE FLIGHT PHASE
...
QUALITY OF BALANCING
TOTALS HorvitzThompson_estimators Relative_deviation
pik 1.000000e+02 1.000000e+02 1.406875e-12

2 1.000000e+03 1.000000e+03 0.000000e+00

3 4.940807e+02 4.942810e+02 4.054567e-02

4 4.948565e+02 4.943465e+02 -1.030694e-01

5 3.301210e+02 3.305846e+02 1.404206e-01

6 3.247878e+02 3.258948e+02 3.408236e-01

7 3.197442e-14 -1.792329e+01 -5.605509e+16

8 6.750156e-14 3.396346e+00 5.031507e+15

9 -8.082424e-14 -3.028730e+01 3.747304e+16