R Graphical Manual

Browse All

Last data update: 2014.03.03

R: A kernel average smoother function to weigh residuals...

gaussianWeights

R Documentation

A kernel average smoother function to weigh residuals according to a Gaussian density function This function is still experimental... use with care

Description

A calibration dataset in database format (cf. modCost for the database format) is extended in order to fit model output using a weighted least squares approach. To this end, the observations are replicated for a certain number of times, and weights are assigned to the replicates according to a Gaussian density function. This density has the relevant observation as mean value. The standard deviation, provided as a parameter, determines the number of inserted replicate observations (see Detail).

This weighted regression approach may be interesting when discontinuities exist in the observational data. Under these circumstances small changes in the timing (or more general the position along the axis of the independent variable) of the model output may have a disproportional impact on the overall goodness-of-fit (e.g. timing of nutrient depletion). Additionally, this approach may be used to model uncertainty in the independent variable (e.g. slices of sediment profiles, or the timing of a sampling).

Usage

gaussianWeights (obs, x = x, y = y, xmodel, spread, weight = "none",
                 aggregation = x ,ordering)

Arguments

`obs`	dataset in long (database) format as is typically used by modCost
`x`	name of the independent variable (typically x, cf. modCost) in `obs`. Defaults to x (not given as character string; cf. subset)
`y`	name of the dependent variable in `obs`. Defaults to y.
`xmodel`	an ordered vector of unique times at which model output is produced. If not given, the independent variable of the observational dataset is used.
`spread`	standard deviation used to calculate the weights from a normal density function. This value also determines the number of points from the model output that are compared to a specific observa- tion in `obs` (2 * 3 * spread + 1; containing 99.7% of the Gaussian distribution, centered around the observation of interest).
`weight`	scaling factor of the modCost function ("sd", "mean", or "none"). The Gaussian weights are multiplied by this factor to account for differences in units.
`aggregation`	vector of column names from the dataset that are used to aggregate observations while calculating the scaling factor. Defaults to the variable name, "name".
`ordering`	Optional extra grouping and ordering of observations. Given as a vector of variable names. If none given, ordering will be done by variable name and independent variable. If both aggregation and ordering variables are given, ordering will be done as follows: x within ordering (in reverse order) within aggregation (in reverse order). Aggregation and ordering should be disjoint sets of variable names.

Details

Suppose: spread = 1/24 (days; = 1 hour) x = time in days, 1 per hour

Then: obs_i is replicated 7 times (spread = observational periodicity = 1 hour):

=> obs_i-3 = ... = obs_i-1 = obs_i = obs_i+1 = ... = obs_i+3

The weights (W_i+j, for j = -3 ...3) are calculated as follows: W'_i+j = 1/(spread * sqrt(2pi)) * exp(-1/2 * ((obs_i+j - obs_i)/spread)^2

W_i+j = W'_i+j/sum(W_i-3,...,W_i+3) (such that their sum equals 1)

Value

A modified version of obs is returned with the following extensions:

1. Each observation obs[i] is replicated n times were n represents the number of modelx values within the interval [obs_i - (3 * spread), obs_i + 3 * spread)].

2. These replicate observations get the same x values as their modeled counterparts (xmodel).

3. Weights are given in column, called "err"

The returned data frame has the following columns:

"name" or another name specified by the first element of aggregation. Usually this column contains the names of the observed variables.
"x" or another name specified by x
"y" or another name specified by y
"err" containing the calculated weights
The rest of the columns of the data frame given by obs in that order.

Author(s)

Tom Van Engeland <tom.vanengeland@nioz.nl>

Examples

## =======================================================================
## A Sediment example
## =======================================================================

## Sediment oxygen concentration is measured every
## centimeter in 3 sediment types
depth <- 0:7
observations <- data.frame(
                    profile = rep(c("mud","silt","sand"), each=8),
                    depth   = depth,
                    O2      = c(c(6,1,0.5,0.1,0.05,0,0,0),
                                c(6,5,3,2,1.5,1,0.5,0),
                                c(6,6,5,4,3,2,1,0)
                              )
                )

## A model generates profiles with a depth resolution of 1 millimeter
modeldepths <- seq(0, 9, by = 0.05)

## All these model outputs are compared with  weighed observations.
gaussianWeights(obs = observations, x = depth, y = O2,
                xmodel = modeldepths,
                spread = 0.1, weight = "none", 
                aggregation = profile)



# Weights of one observation in silt at depth 2:
Sub <- subset(observations, subset = (profile == "silt" & depth == 2))
plot(Sub[,-1])
SubWW <- gaussianWeights(obs = Sub, x = depth, y = O2, 
                xmodel = modeldepths, spread = 0.5, 
                weight="none", aggregation = profile)
SubWW[,-1]

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(FME)
Loading required package: deSolve

Attaching package: 'deSolve'

The following object is masked from 'package:graphics':

    matplot

Loading required package: rootSolve
Loading required package: coda
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/FME/gaussianWeights.Rd_%03d_medium.png", width=480, height=480)
> ### Name: gaussianWeights
> ### Title: A kernel average smoother function to weigh residuals according
> ###   to a Gaussian density function This function is still experimental...
> ###   use with care
> ### Aliases: gaussianWeights
> ### Keywords: utilities
> 
> ### ** Examples
> 
> ## =======================================================================
> ## A Sediment example
> ## =======================================================================
> 
> ## Sediment oxygen concentration is measured every
> ## centimeter in 3 sediment types
> depth <- 0:7
> observations <- data.frame(
+                     profile = rep(c("mud","silt","sand"), each=8),
+                     depth   = depth,
+                     O2      = c(c(6,1,0.5,0.1,0.05,0,0,0),
+                                 c(6,5,3,2,1.5,1,0.5,0),
+                                 c(6,6,5,4,3,2,1,0)
+                               )
+                 )
> 
> ## A model generates profiles with a depth resolution of 1 millimeter
> modeldepths <- seq(0, 9, by = 0.05)
> 
> ## All these model outputs are compared with  weighed observations.
> gaussianWeights(obs = observations, x = depth, y = O2,
+                 xmodel = modeldepths,
+                 spread = 0.1, weight = "none", 
+                 aggregation = profile)
    profile depth   O2        err
1       mud  0.00 6.00   3.004061
2       mud  0.05 6.00   3.404047
3       mud  0.10 6.00   4.952860
4       mud  0.15 6.00   9.253160
5       mud  0.20 6.00  22.197177
6       mud  0.25 6.00  68.372119
7       mud  0.30 6.00 270.416975
22      mud  0.70 1.00 450.816819
23      mud  0.75 1.00 113.984342
24      mud  0.80 1.00  37.005298
25      mud  0.85 1.00  15.426103
26      mud  0.90 1.00   8.256998
27      mud  0.95 1.00   5.674946
28      mud  1.00 1.00   5.008122
29      mud  1.05 1.00   5.674946
30      mud  1.10 1.00   8.256998
31      mud  1.15 1.00  15.426103
32      mud  1.20 1.00  37.005298
33      mud  1.25 1.00 113.984342
34      mud  1.30 1.00 450.816819
61      mud  1.70 0.50 449.816819
62      mud  1.75 0.50 113.731503
63      mud  1.80 0.50  36.923213
64      mud  1.85 0.50  15.391885
65      mud  1.90 0.50   8.238682
66      mud  1.95 0.50   5.662358
67      mud  2.00 0.50   4.997013
68      mud  2.05 0.50   5.662358
69      mud  2.10 0.50   8.238682
70      mud  2.15 0.50  15.391885
71      mud  2.20 0.50  36.923213
72      mud  2.25 0.50 113.731503
97      mud  2.70 0.10 449.816819
98      mud  2.75 0.10 113.731503
99      mud  2.80 0.10  36.923213
100     mud  2.85 0.10  15.391885
101     mud  2.90 0.10   8.238682
102     mud  2.95 0.10   5.662358
103     mud  3.00 0.10   4.997013
104     mud  3.05 0.10   5.662358
105     mud  3.10 0.10   8.238682
106     mud  3.15 0.10  15.391885
107     mud  3.20 0.10  36.923213
108     mud  3.25 0.10 113.731503
133     mud  3.70 0.05 450.816819
134     mud  3.75 0.05 113.984342
135     mud  3.80 0.05  37.005298
136     mud  3.85 0.05  15.426103
137     mud  3.90 0.05   8.256998
138     mud  3.95 0.05   5.674946
139     mud  4.00 0.05   5.008122
140     mud  4.05 0.05   5.674946
141     mud  4.10 0.05   8.256998
142     mud  4.15 0.05  15.426103
143     mud  4.20 0.05  37.005298
144     mud  4.25 0.05 113.984342
145     mud  4.30 0.05 450.816819
172     mud  4.70 0.00 449.816819
173     mud  4.75 0.00 113.731503
174     mud  4.80 0.00  36.923213
175     mud  4.85 0.00  15.391885
176     mud  4.90 0.00   8.238682
177     mud  4.95 0.00   5.662358
178     mud  5.00 0.00   4.997013
179     mud  5.05 0.00   5.662358
180     mud  5.10 0.00   8.238682
181     mud  5.15 0.00  15.391885
182     mud  5.20 0.00  36.923213
183     mud  5.25 0.00 113.731503
208     mud  5.70 0.00 449.816819
209     mud  5.75 0.00 113.731503
210     mud  5.80 0.00  36.923213
211     mud  5.85 0.00  15.391885
212     mud  5.90 0.00   8.238682
213     mud  5.95 0.00   5.662358
214     mud  6.00 0.00   4.997013
215     mud  6.05 0.00   5.662358
216     mud  6.10 0.00   8.238682
217     mud  6.15 0.00  15.391885
218     mud  6.20 0.00  36.923213
219     mud  6.25 0.00 113.731503
244     mud  6.70 0.00 449.816819
245     mud  6.75 0.00 113.731503
246     mud  6.80 0.00  36.923213
247     mud  6.85 0.00  15.391885
248     mud  6.90 0.00   8.238682
249     mud  6.95 0.00   5.662358
250     mud  7.00 0.00   4.997013
251     mud  7.05 0.00   5.662358
252     mud  7.10 0.00   8.238682
253     mud  7.15 0.00  15.391885
254     mud  7.20 0.00  36.923213
255     mud  7.25 0.00 113.731503
15     sand  0.00 6.00   3.004061
16     sand  0.05 6.00   3.404047
17     sand  0.10 6.00   4.952860
18     sand  0.15 6.00   9.253160
19     sand  0.20 6.00  22.197177
20     sand  0.25 6.00  68.372119
21     sand  0.30 6.00 270.416975
48     sand  0.70 6.00 450.816819
49     sand  0.75 6.00 113.984342
50     sand  0.80 6.00  37.005298
51     sand  0.85 6.00  15.426103
52     sand  0.90 6.00   8.256998
53     sand  0.95 6.00   5.674946
54     sand  1.00 6.00   5.008122
55     sand  1.05 6.00   5.674946
56     sand  1.10 6.00   8.256998
57     sand  1.15 6.00  15.426103
58     sand  1.20 6.00  37.005298
59     sand  1.25 6.00 113.984342
60     sand  1.30 6.00 450.816819
85     sand  1.70 5.00 449.816819
86     sand  1.75 5.00 113.731503
87     sand  1.80 5.00  36.923213
88     sand  1.85 5.00  15.391885
89     sand  1.90 5.00   8.238682
90     sand  1.95 5.00   5.662358
91     sand  2.00 5.00   4.997013
92     sand  2.05 5.00   5.662358
93     sand  2.10 5.00   8.238682
94     sand  2.15 5.00  15.391885
95     sand  2.20 5.00  36.923213
96     sand  2.25 5.00 113.731503
121    sand  2.70 4.00 449.816819
122    sand  2.75 4.00 113.731503
123    sand  2.80 4.00  36.923213
124    sand  2.85 4.00  15.391885
125    sand  2.90 4.00   8.238682
126    sand  2.95 4.00   5.662358
127    sand  3.00 4.00   4.997013
128    sand  3.05 4.00   5.662358
129    sand  3.10 4.00   8.238682
130    sand  3.15 4.00  15.391885
131    sand  3.20 4.00  36.923213
132    sand  3.25 4.00 113.731503
159    sand  3.70 3.00 450.816819
160    sand  3.75 3.00 113.984342
161    sand  3.80 3.00  37.005298
162    sand  3.85 3.00  15.426103
163    sand  3.90 3.00   8.256998
164    sand  3.95 3.00   5.674946
165    sand  4.00 3.00   5.008122
166    sand  4.05 3.00   5.674946
167    sand  4.10 3.00   8.256998
168    sand  4.15 3.00  15.426103
169    sand  4.20 3.00  37.005298
170    sand  4.25 3.00 113.984342
171    sand  4.30 3.00 450.816819
196    sand  4.70 2.00 449.816819
197    sand  4.75 2.00 113.731503
198    sand  4.80 2.00  36.923213
199    sand  4.85 2.00  15.391885
200    sand  4.90 2.00   8.238682
201    sand  4.95 2.00   5.662358
202    sand  5.00 2.00   4.997013
203    sand  5.05 2.00   5.662358
204    sand  5.10 2.00   8.238682
205    sand  5.15 2.00  15.391885
206    sand  5.20 2.00  36.923213
207    sand  5.25 2.00 113.731503
232    sand  5.70 1.00 449.816819
233    sand  5.75 1.00 113.731503
234    sand  5.80 1.00  36.923213
235    sand  5.85 1.00  15.391885
236    sand  5.90 1.00   8.238682
237    sand  5.95 1.00   5.662358
238    sand  6.00 1.00   4.997013
239    sand  6.05 1.00   5.662358
240    sand  6.10 1.00   8.238682
241    sand  6.15 1.00  15.391885
242    sand  6.20 1.00  36.923213
243    sand  6.25 1.00 113.731503
268    sand  6.70 0.00 449.816819
269    sand  6.75 0.00 113.731503
270    sand  6.80 0.00  36.923213
271    sand  6.85 0.00  15.391885
272    sand  6.90 0.00   8.238682
273    sand  6.95 0.00   5.662358
274    sand  7.00 0.00   4.997013
275    sand  7.05 0.00   5.662358
276    sand  7.10 0.00   8.238682
277    sand  7.15 0.00  15.391885
278    sand  7.20 0.00  36.923213
279    sand  7.25 0.00 113.731503
8      silt  0.00 6.00   3.004061
9      silt  0.05 6.00   3.404047
10     silt  0.10 6.00   4.952860
11     silt  0.15 6.00   9.253160
12     silt  0.20 6.00  22.197177
13     silt  0.25 6.00  68.372119
14     silt  0.30 6.00 270.416975
35     silt  0.70 5.00 450.816819
36     silt  0.75 5.00 113.984342
37     silt  0.80 5.00  37.005298
38     silt  0.85 5.00  15.426103
39     silt  0.90 5.00   8.256998
40     silt  0.95 5.00   5.674946
41     silt  1.00 5.00   5.008122
42     silt  1.05 5.00   5.674946
43     silt  1.10 5.00   8.256998
44     silt  1.15 5.00  15.426103
45     silt  1.20 5.00  37.005298
46     silt  1.25 5.00 113.984342
47     silt  1.30 5.00 450.816819
73     silt  1.70 3.00 449.816819
74     silt  1.75 3.00 113.731503
75     silt  1.80 3.00  36.923213
76     silt  1.85 3.00  15.391885
77     silt  1.90 3.00   8.238682
78     silt  1.95 3.00   5.662358
79     silt  2.00 3.00   4.997013
80     silt  2.05 3.00   5.662358
81     silt  2.10 3.00   8.238682
82     silt  2.15 3.00  15.391885
83     silt  2.20 3.00  36.923213
84     silt  2.25 3.00 113.731503
109    silt  2.70 2.00 449.816819
110    silt  2.75 2.00 113.731503
111    silt  2.80 2.00  36.923213
112    silt  2.85 2.00  15.391885
113    silt  2.90 2.00   8.238682
114    silt  2.95 2.00   5.662358
115    silt  3.00 2.00   4.997013
116    silt  3.05 2.00   5.662358
117    silt  3.10 2.00   8.238682
118    silt  3.15 2.00  15.391885
119    silt  3.20 2.00  36.923213
120    silt  3.25 2.00 113.731503
146    silt  3.70 1.50 450.816819
147    silt  3.75 1.50 113.984342
148    silt  3.80 1.50  37.005298
149    silt  3.85 1.50  15.426103
150    silt  3.90 1.50   8.256998
151    silt  3.95 1.50   5.674946
152    silt  4.00 1.50   5.008122
153    silt  4.05 1.50   5.674946
154    silt  4.10 1.50   8.256998
155    silt  4.15 1.50  15.426103
156    silt  4.20 1.50  37.005298
157    silt  4.25 1.50 113.984342
158    silt  4.30 1.50 450.816819
184    silt  4.70 1.00 449.816819
185    silt  4.75 1.00 113.731503
186    silt  4.80 1.00  36.923213
187    silt  4.85 1.00  15.391885
188    silt  4.90 1.00   8.238682
189    silt  4.95 1.00   5.662358
190    silt  5.00 1.00   4.997013
191    silt  5.05 1.00   5.662358
192    silt  5.10 1.00   8.238682
193    silt  5.15 1.00  15.391885
194    silt  5.20 1.00  36.923213
195    silt  5.25 1.00 113.731503
220    silt  5.70 0.50 449.816819
221    silt  5.75 0.50 113.731503
222    silt  5.80 0.50  36.923213
223    silt  5.85 0.50  15.391885
224    silt  5.90 0.50   8.238682
225    silt  5.95 0.50   5.662358
226    silt  6.00 0.50   4.997013
227    silt  6.05 0.50   5.662358
228    silt  6.10 0.50   8.238682
229    silt  6.15 0.50  15.391885
230    silt  6.20 0.50  36.923213
231    silt  6.25 0.50 113.731503
256    silt  6.70 0.00 449.816819
257    silt  6.75 0.00 113.731503
258    silt  6.80 0.00  36.923213
259    silt  6.85 0.00  15.391885
260    silt  6.90 0.00   8.238682
261    silt  6.95 0.00   5.662358
262    silt  7.00 0.00   4.997013
263    silt  7.05 0.00   5.662358
264    silt  7.10 0.00   8.238682
265    silt  7.15 0.00  15.391885
266    silt  7.20 0.00  36.923213
267    silt  7.25 0.00 113.731503
> 
> 
> 
> # Weights of one observation in silt at depth 2:
> Sub <- subset(observations, subset = (profile == "silt" & depth == 2))
> plot(Sub[,-1])
> SubWW <- gaussianWeights(obs = Sub, x = depth, y = O2, 
+                 xmodel = modeldepths, spread = 0.5, 
+                 weight="none", aggregation = profile)
> SubWW[,-1]
   depth O2        err
1   0.50  3 2251.25311
2   0.55  3 1676.12905
3   0.60  3 1260.47294
4   0.65  3  957.42011
5   0.70  3  734.53840
6   0.75  3  569.20593
7   0.80  3  445.52006
8   0.85  3  352.21517
9   0.90  3  281.24953
10  0.95  3  226.83940
11  1.00  3  184.79411
12  1.05  3  152.05500
13  1.10  3  126.37356
14  1.15  3  106.08517
15  1.20  3   89.94895
16  1.25  3   77.03365
17  1.30  3   66.63582
18  1.35  3   58.22078
19  1.40  3   51.37966
20  1.45  3   45.79808
21  1.50  3   41.23314
22  1.55  3   37.49630
23  1.60  3   34.44081
24  1.65  3   31.95224
25  1.70  3   29.94140
26  1.75  3   28.33909
27  1.80  3   27.09210
28  1.85  3   26.16028
29  1.90  3   25.51438
30  1.95  3   25.13452
31  2.00  3   25.00916
32  2.05  3   25.13452
33  2.10  3   25.51438
34  2.15  3   26.16028
35  2.20  3   27.09210
36  2.25  3   28.33909
37  2.30  3   29.94140
38  2.35  3   31.95224
39  2.40  3   34.44081
40  2.45  3   37.49630
41  2.50  3   41.23314
42  2.55  3   45.79808
43  2.60  3   51.37966
44  2.65  3   58.22078
45  2.70  3   66.63582
46  2.75  3   77.03365
47  2.80  3   89.94895
48  2.85  3  106.08517
49  2.90  3  126.37356
50  2.95  3  152.05500
51  3.00  3  184.79411
52  3.05  3  226.83940
53  3.10  3  281.24953
54  3.15  3  352.21517
55  3.20  3  445.52006
56  3.25  3  569.20593
57  3.30  3  734.53840
58  3.35  3  957.42011
59  3.40  3 1260.47294
60  3.45  3 1676.12905
61  3.50  3 2251.25311
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>