/
cov_panel_dsgn.R
157 lines (146 loc) · 6.21 KB
/
cov_panel_dsgn.R
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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
###############################################################################
# Function: cov_panel_dsgn (exported)
# Programmers: Tom Kincaid
# Tony Olsen
# Date: March 14, 2019
#
#'
#' Create a covariance matrix for a panel design
#'
#' Covariance structure accounts for the panel design and the four variance
#' components: unit variation, period variation, unit by period interaction
#' variation and index (or residual) variation. The model incorporates unit,
#' period, unit by period, and index variance components. It also includes a
#' provision for unit correlation and period autocorrelation.
#'
#' @param paneldsgn A matrix (dimensions: number of panels (rows) by number of
#' periods (columns)) containing the number of units visited for each
#' combination of panel and period. Default is matrix(50, 1, 10) which is a
#' single panel of 50 units visited 10 times, typical time is a period.
#'
#' @param nrepeats Either \code{NULL} or a list of matrices the same length as
#' paneldsgn specifying the number of revisits made to units in a panel in the
#' same period for each design. Specifying \code{NULL} indicates that number of
#' revisits to units is the same for all panels and for all periods and for
#' all panel designs. The default is \code{NULL}, a single visit. Names must match
#' list names in \code{paneldsgn}.
#'
#' @param unit_var The variance component estimate for unit. The default is
#' \code{NULL}.
#'
#' @param period_var The variance component estimate for period The default is
#' \code{NULL}.
#'
#' @param unitperiod_var The variance component estimate for unit by period
#' interaction. The default is \code{NULL}.
#'
#' @param index_var The variance component estimate for index error. The
#' default is \code{NULL}.
#'
#' @param unit_rho Unit correlation across periods. The default is \code{1}.
#'
#' @param period_rho Period autocorrelation. The default is \code{0}.
#'
#' @details Covariance structure accounts for the panel design and the four
#' variance components: unit variation, period variation, unit by period
#' interaction variation and index (or residual) variation. Uses the model
#' structure defined by Urquhart 2012.
#'
#' If \code{nrepeats} is \code{NULL}, then no units sampled more than once in a specific
#' panel, period combination) and then unit by period and index variances are
#' added together or user may have only estimated unit, period and unit by
#' period variance components so that index component is zero. It calculates
#' the covariance matrix for the simple linear regression. The standard error
#' for a linear trend coefficient is the square root of the variance.
#'
#' @references
#' Urquhart, N. S., W. S. Overton, et al. (1993) Comparing sampling designs
#' for monitoring ecological status and trends: impact of temporal patterns.
#' In: \emph{Statistics for the Environment.} V. Barnett and K. F. Turkman.
#' John Wiley & Sons, New York, pp. 71-86.
#'
#' Urquhart, N. S. and T. M. Kincaid (1999). Designs for detecting trends
#' from repeated surveys of ecological resources. \emph{Journal of
#' Agricultural, Biological, and Environmental Statistics}, \bold{4(4)},
#' 404-414.
#'
#' Urquhart, N. S. (2012). The role of monitoring design in detecting trend in
#' long-term ecological monitoring studies. In: \emph{Design and Analysis of
#' Long-term Ecological Monitoring Studies.} R. A. Gitzen, J. J. Millspaugh,
#' A. B. Cooper, and D. S. Licht (eds.). Cambridge University Press, New York,
#' pp. 151-173.
#
#' @return A list containing the covariance matrix (\code{cov}) for the panel design,
#' the input panel design (\code{paneldsgn}), the input \code{nrepeats} design
#' (\code{nrepeats.dsgn}) and the function call.
#'
#' @author Tony Olsen \email{Olsen.Tony@@epa.gov}
#'
#' @seealso
#' \describe{
#' \item{\code{\link{power_dsgn}}}{ for power calculations of multiple panel
#' designs}
#' }
#'
#' @keywords survey
#'
#' @export
###############################################################################
cov_panel_dsgn <- function(paneldsgn = matrix(50, 1, 10), nrepeats = 1,
unit_var = NULL, period_var = NULL, unitperiod_var = NULL, index_var = NULL,
unit_rho = 1, period_rho = 0) {
paneldsgn <- as.matrix(paneldsgn)
nperiod <- ncol(paneldsgn)
npanel <- nrow(paneldsgn)
if (any(is.null(unit_var), is.null(period_var), is.null(unitperiod_var), is.null(index_var))) {
stop("Must provide four variance components. index_var may be zero if nrepeats=1")
}
# Create covariance matrix for period correlation
rhomat <- matrix(0, nrow = nperiod, ncol = nperiod)
if (nperiod > 1) {
rhomat[1, ] <- 0:(nperiod - 1)
for (i in 2:nperiod) {
rhomat[i, ] <- c(rev(1:(i - 1)), 0:(nperiod - i))
}
}
period_cov <- period_var * (period_rho^rhomat)
# Create covariance matrix for unit correlation
unit_cov <- unit_var * (unit_rho^rhomat)
# Construct panel covariance term
# nunits is number of units in a panel
nunits <- apply(paneldsgn, 1, max)
if (npanel > 1) {
pan_cov <- kronecker(unit_cov, diag(ifelse(nunits > 0, 1 / nunits, 0)))
} else {
pan_cov <- unit_cov / nunits
}
# construct period covariance term
yr_cov <- kronecker(period_cov, matrix(1, nrow = npanel, ncol = npanel))
# Construct unit by period covariance term
if (length(nunits) > 1) {
sy_cov <- unitperiod_var * kronecker(diag(nperiod), diag(ifelse(nunits > 0, 1 / nunits, 0)))
} else {
sy_cov <- unitperiod_var * diag(nperiod) / nunits
}
# Construct index covariance term
# Create nrepeats revisit design if necessary
if (is.null(nrepeats)) {
nrepeats <- vector("list", length(paneldsgn))
for (i in names(paneldsgn)) {
nrepeats[[i]] <- paneldsgn[[i]]
nrepeats[[i]][nrepeats[[i]] > 0] <- 1
}
}
if (npanel > 1) {
index_cov <- index_var * kronecker(
diag(nperiod),
diag(ifelse(nunits > 0, 1 / (nunits * apply(nrepeats, 1, max)), 0))
)
} else {
index_cov <- index_var * diag(nperiod) / nunits
}
# overall covariance matrix
phi <- pan_cov + yr_cov + sy_cov + index_cov
cov <- list(cov = phi, paneldsgn = paneldsgn, nrepeat_dsgn = nrepeats, call = sys.call())
return(cov)
}