Stem.Bootstrap: Parametric bootstrap
In Stem: Spatio-temporal models in R
Description
Usage
Arguments
Details
Value
Author(s)
References
See Also
Examples
Description
This functions performs the spatio-temporal parametric bootstrap for computing the parameter standard errors.
Usage
1
Stem.Bootstrap(StemModel, B)
Arguments
StemModel
an object of class “Stem.Model” given as output by the Stem.Estimation function.
B
number of bootstrap iterations.
Details
The spatio-temporal bootstrap is used for parameter uncertainty assessment. The resampling scheme is based on the estimated model: each bootstrap sample is drawn
directly from the Gaussian distributions which define the model with the parameter vector replaced by the corresponding ML estimates. For each of the B bootstrap
samples, the ML estimates are computed (using the procedure of Stem.Estimation function). Then the B bootstrap replications are returned in a list.
Value
The function returns a list of elements called “boot.output”. Each element of the list is an object of class “Stem.Model” and so
it is composed by the following elements:
skeleton
a list with components phi, p, K
where the phi vector is given by
the StemModel$estimates$phi.hat vector.
data
a list with components z, coordinates, covariates, r, n and d where the z matrix is
given by the simulated data matrix.
estimates
A list of four objects: phi.hat, y.smoothed, loglik, convergence.par as the output of the Stem.Estimation function.
Author(s)
Michela Cameletti michela.cameletti@unibg.it
References
Amisigo, B.A., Van De Giesen, N.C. (2005) Using a spatio-temporal dynamic state-space model with the EM algorithm to patch gaps in daily riverflow series. Hydrology and Earth System Sciences 9, 209–224.
Fasso', A., Cameletti, M., Nicolis, O. (2007) Air quality monitoring using heterogeneous networks. Environmetrics 18, 245–264.
Fasso', A., Cameletti, M. (2007) A general spatio-temporal model for environmental data. Tech.rep. n.27 Graspa - The Italian Group of Environmental Statistics - http://www.graspa.org.
Fasso', A., Cameletti, M. (2009) A unified statistical approach for simulation, modelling, analysis and mapping of environmental data. Accepted for publication by Simulation: transaction of the Society for Modeling and Simulation International.
Xu, K., Wikle, C.K. (2007) Estimation of parameterized spatio-temporal dynamic models. Journal of Statistical Inference and Planning 137, 567–588.
Mc Lachlan, G.J., Krishnan, T. (1997) The EM Algorithm and Extensions. Wiley, New York.
Shumway, R.H., Stoffer, D.S. (2006) Time Series Analysis and Its Applications: with R Examples. Springer, New York.
See Also
See Also pm10, Stem.Model and Stem.Estimation
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#load the data
data(pm10)
#extract the data
coordinates <- pm10$coords
covariates <- pm10$covariates
z <- pm10$z
#build the parameter list
#(the phi list is used for the algorithm starting values)
phi <- list(beta=matrix(c(3.65,0.046,-0.904),3,1),
sigma2eps=0.1,
sigma2omega=0.2,
theta=0.01,
G=matrix(0.77,1,1),
Sigmaeta=matrix(0.3,1,1),
m0=as.matrix(0),
C0=as.matrix(1))
K <- matrix(1,ncol(z),1)
mod1 <- Stem.Model(z=z,covariates=covariates,
coordinates=coordinates,phi=phi,K=K)
class(mod1)
is.Stem.Model(mod1)
#mod1 is given as output by the Stem.Model function
mod1.est <- Stem.Estimation(mod1)
#it is computer intensive
mod1.boot <- Stem.Bootstrap(StemModel=mod1.est, B=3)
names(mod1.boot)
#the first element of the output list
names(mod1.boot$boot.output[[1]])
#check if there is no convergence for some bootstrap iteration
B <- length(mod1.boot$boot.output)
n.null <- sum (unlist (lapply(mod1.boot$boot.output,
function(x) length(x$StemModel)==1)) )
pos.null <- which((unlist (lapply(mod1.boot$boot.output,
function(x) length(x$StemModel)==1)) ))
cat("------------B non null =", B - n.null,"\n")
if(length(pos.null)>0) boot1.mod = mod1.boot[-pos.null]
#put the bootstrap output in a matrix
npar <- length(unlist((mod1.boot$boot.output[[1]]$estimates$phi.hat)))-1
boot.estimates <- matrix(NA, nrow = (B - n.null), ncol = npar)
for(b in 1:(B - n.null)) {
phi.estimated <- mod1.boot$boot.output[[b]]$estimates$phi.hat
boot.estimates[b,] <- c(phi.estimated$beta,
phi.estimated$sigma2eps,
phi.estimated$sigma2omega,
phi.estimated$theta,
phi.estimated$G,
phi.estimated$Sigmaeta,
phi.estimated$m0)
}
#compute the parameter standard errors
se <- sqrt(diag(var(na.omit(boot.estimates))))
#create a summary table with Estimates, Standard Errors (SE) and T-statistics.
phi.hat <- mod1.est$estimates$phi.hat
MLE <- c(phi.hat$beta, phi.hat$sigma2eps, phi.hat$sigma2omega,
phi.hat$theta, phi.hat$G, phi.hat$Sigmaeta,phi.hat$m0)
output1 <- cbind(MLE, se, MLE/se)
colnames(output1)<- c("Estimate", "SE", "T-stat.")
output1
#compute the 95% confidence intervals
IC <- matrix(NA,nrow=npar,ncol=2)
for(i in 1 : npar) {
IC[i,] <- c(quantile(boot.estimates[,i],0.025),
quantile(boot.estimates[,i],0.975))
}
#create a summary table with Estimates, Standard Errors (SE)
#and T-statistics and confidence intervals.
output2 <- cbind(output1,IC)
colnames(output2) <- c("Estimate", "SE", "T-stat.", "IC_inf", "IC_sup")
output2
Example output
Loading required package: mvtnorm
Loading required package: MASS
[1] "Stem.Model"
[1] TRUE
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
***NR Algorithm - iteration n. 3
***NR Algorithm - iteration n. 4
***NR Algorithm - iteration n. 5
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
***NR Algorithm - iteration n. 3
***NR Algorithm - iteration n. 4
****************EM Algorithm - iteration n. 3
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
***NR Algorithm - iteration n. 3
***NR Algorithm - iteration n. 4
***NR Algorithm - iteration n. 5
****************EM Algorithm - iteration n. 4
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 5
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 6
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 7
=========================== BOOTSTRAP ITERATION N. 1 ===========================
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 3
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 4
=========================== BOOTSTRAP ITERATION N. 2 ===========================
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 3
=========================== BOOTSTRAP ITERATION N. 3 ===========================
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 3
[1] "boot.output"
[1] "skeleton" "data" "estimates"
------------B non null = 3
Estimate SE T-stat.
[1,] 4.783145941 1.562749e-02 306.072608
[2,] 0.002709126 8.723504e-05 31.055474
[3,] -0.947914091 3.846030e-02 -24.646557
[4,] 0.084248655 2.542421e-03 33.137183
[5,] 0.173456019 1.713464e-03 101.231229
[6,] 0.007540571 3.808402e-04 19.799832
[7,] 0.975095649 7.088034e-03 137.569271
[8,] 0.067555184 1.922474e-03 35.139718
[9,] -1.255130285 4.492622e-01 -2.793759
Estimate SE T-stat. IC_inf IC_sup
[1,] 4.783145941 1.562749e-02 306.072608 4.777880916 4.805210495
[2,] 0.002709126 8.723504e-05 31.055474 0.002715966 0.002870147
[3,] -0.947914091 3.846030e-02 -24.646557 -1.003994089 -0.940250813
[4,] 0.084248655 2.542421e-03 33.137183 0.081054309 0.085879295
[5,] 0.173456019 1.713464e-03 101.231229 0.174186044 0.177415693
[6,] 0.007540571 3.808402e-04 19.799832 0.007524053 0.008226455
[7,] 0.975095649 7.088034e-03 137.569271 0.972485276 0.984340910
[8,] 0.067555184 1.922474e-03 35.139718 0.068947486 0.072599600
[9,] -1.255130285 4.492622e-01 -2.793759 -0.922237838 -0.072558155
Stem documentation built on May 2, 2019, 8:56 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Author(s) References See Also Examples
Description
This functions performs the spatio-temporal parametric bootstrap for computing the parameter standard errors.
Usage
1 | Stem.Bootstrap(StemModel, B)
|
Arguments
StemModel |
an object of class “Stem.Model” given as output by the |
B |
number of bootstrap iterations. |
Details
The spatio-temporal bootstrap is used for parameter uncertainty assessment. The resampling scheme is based on the estimated model: each bootstrap sample is drawn
directly from the Gaussian distributions which define the model with the parameter vector replaced by the corresponding ML estimates. For each of the B bootstrap
samples, the ML estimates are computed (using the procedure of Stem.Estimation function). Then the B bootstrap replications are returned in a list.
Value
The function returns a list of elements called “boot.output”. Each element of the list is an object of class “Stem.Model” and so it is composed by the following elements:
skeleton |
a list with components |
data |
a list with components |
estimates |
A list of four objects: |
Author(s)
Michela Cameletti michela.cameletti@unibg.it
References
Amisigo, B.A., Van De Giesen, N.C. (2005) Using a spatio-temporal dynamic state-space model with the EM algorithm to patch gaps in daily riverflow series. Hydrology and Earth System Sciences 9, 209–224.
Fasso', A., Cameletti, M., Nicolis, O. (2007) Air quality monitoring using heterogeneous networks. Environmetrics 18, 245–264.
Fasso', A., Cameletti, M. (2007) A general spatio-temporal model for environmental data. Tech.rep. n.27 Graspa - The Italian Group of Environmental Statistics - http://www.graspa.org.
Fasso', A., Cameletti, M. (2009) A unified statistical approach for simulation, modelling, analysis and mapping of environmental data. Accepted for publication by Simulation: transaction of the Society for Modeling and Simulation International.
Xu, K., Wikle, C.K. (2007) Estimation of parameterized spatio-temporal dynamic models. Journal of Statistical Inference and Planning 137, 567–588.
Mc Lachlan, G.J., Krishnan, T. (1997) The EM Algorithm and Extensions. Wiley, New York.
Shumway, R.H., Stoffer, D.S. (2006) Time Series Analysis and Its Applications: with R Examples. Springer, New York.
See Also
See Also pm10, Stem.Model and Stem.Estimation
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 | #load the data
data(pm10)
#extract the data
coordinates <- pm10$coords
covariates <- pm10$covariates
z <- pm10$z
#build the parameter list
#(the phi list is used for the algorithm starting values)
phi <- list(beta=matrix(c(3.65,0.046,-0.904),3,1),
sigma2eps=0.1,
sigma2omega=0.2,
theta=0.01,
G=matrix(0.77,1,1),
Sigmaeta=matrix(0.3,1,1),
m0=as.matrix(0),
C0=as.matrix(1))
K <- matrix(1,ncol(z),1)
mod1 <- Stem.Model(z=z,covariates=covariates,
coordinates=coordinates,phi=phi,K=K)
class(mod1)
is.Stem.Model(mod1)
#mod1 is given as output by the Stem.Model function
mod1.est <- Stem.Estimation(mod1)
#it is computer intensive
mod1.boot <- Stem.Bootstrap(StemModel=mod1.est, B=3)
names(mod1.boot)
#the first element of the output list
names(mod1.boot$boot.output[[1]])
#check if there is no convergence for some bootstrap iteration
B <- length(mod1.boot$boot.output)
n.null <- sum (unlist (lapply(mod1.boot$boot.output,
function(x) length(x$StemModel)==1)) )
pos.null <- which((unlist (lapply(mod1.boot$boot.output,
function(x) length(x$StemModel)==1)) ))
cat("------------B non null =", B - n.null,"\n")
if(length(pos.null)>0) boot1.mod = mod1.boot[-pos.null]
#put the bootstrap output in a matrix
npar <- length(unlist((mod1.boot$boot.output[[1]]$estimates$phi.hat)))-1
boot.estimates <- matrix(NA, nrow = (B - n.null), ncol = npar)
for(b in 1:(B - n.null)) {
phi.estimated <- mod1.boot$boot.output[[b]]$estimates$phi.hat
boot.estimates[b,] <- c(phi.estimated$beta,
phi.estimated$sigma2eps,
phi.estimated$sigma2omega,
phi.estimated$theta,
phi.estimated$G,
phi.estimated$Sigmaeta,
phi.estimated$m0)
}
#compute the parameter standard errors
se <- sqrt(diag(var(na.omit(boot.estimates))))
#create a summary table with Estimates, Standard Errors (SE) and T-statistics.
phi.hat <- mod1.est$estimates$phi.hat
MLE <- c(phi.hat$beta, phi.hat$sigma2eps, phi.hat$sigma2omega,
phi.hat$theta, phi.hat$G, phi.hat$Sigmaeta,phi.hat$m0)
output1 <- cbind(MLE, se, MLE/se)
colnames(output1)<- c("Estimate", "SE", "T-stat.")
output1
#compute the 95% confidence intervals
IC <- matrix(NA,nrow=npar,ncol=2)
for(i in 1 : npar) {
IC[i,] <- c(quantile(boot.estimates[,i],0.025),
quantile(boot.estimates[,i],0.975))
}
#create a summary table with Estimates, Standard Errors (SE)
#and T-statistics and confidence intervals.
output2 <- cbind(output1,IC)
colnames(output2) <- c("Estimate", "SE", "T-stat.", "IC_inf", "IC_sup")
output2
|
Example output
Loading required package: mvtnorm
Loading required package: MASS
[1] "Stem.Model"
[1] TRUE
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
***NR Algorithm - iteration n. 3
***NR Algorithm - iteration n. 4
***NR Algorithm - iteration n. 5
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
***NR Algorithm - iteration n. 3
***NR Algorithm - iteration n. 4
****************EM Algorithm - iteration n. 3
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
***NR Algorithm - iteration n. 3
***NR Algorithm - iteration n. 4
***NR Algorithm - iteration n. 5
****************EM Algorithm - iteration n. 4
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 5
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 6
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 7
=========================== BOOTSTRAP ITERATION N. 1 ===========================
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
***NR Algorithm - iteration n. 2
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 3
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 4
=========================== BOOTSTRAP ITERATION N. 2 ===========================
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 3
=========================== BOOTSTRAP ITERATION N. 3 ===========================
****************EM Algorithm - iteration n. 1
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 2
***NR Algorithm - iteration n. 1
****************EM Algorithm - iteration n. 3
[1] "boot.output"
[1] "skeleton" "data" "estimates"
------------B non null = 3
Estimate SE T-stat.
[1,] 4.783145941 1.562749e-02 306.072608
[2,] 0.002709126 8.723504e-05 31.055474
[3,] -0.947914091 3.846030e-02 -24.646557
[4,] 0.084248655 2.542421e-03 33.137183
[5,] 0.173456019 1.713464e-03 101.231229
[6,] 0.007540571 3.808402e-04 19.799832
[7,] 0.975095649 7.088034e-03 137.569271
[8,] 0.067555184 1.922474e-03 35.139718
[9,] -1.255130285 4.492622e-01 -2.793759
Estimate SE T-stat. IC_inf IC_sup
[1,] 4.783145941 1.562749e-02 306.072608 4.777880916 4.805210495
[2,] 0.002709126 8.723504e-05 31.055474 0.002715966 0.002870147
[3,] -0.947914091 3.846030e-02 -24.646557 -1.003994089 -0.940250813
[4,] 0.084248655 2.542421e-03 33.137183 0.081054309 0.085879295
[5,] 0.173456019 1.713464e-03 101.231229 0.174186044 0.177415693
[6,] 0.007540571 3.808402e-04 19.799832 0.007524053 0.008226455
[7,] 0.975095649 7.088034e-03 137.569271 0.972485276 0.984340910
[8,] 0.067555184 1.922474e-03 35.139718 0.068947486 0.072599600
[9,] -1.255130285 4.492622e-01 -2.793759 -0.922237838 -0.072558155
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.