Description Usage Arguments Details Value Author(s) See Also Examples
Produces a bivariate residual plot with simulation polygons to assess goodness-of-fit of bivariate statistical models, provided the user supplies three functions: one to obtain model diagnostics, one to simulate data from a fitted model object, and one to refit the model to simulated data.
1 2 3 4 5 6 7 |
obj |
fitted model object |
sim |
number of simulations used to compute envelope. Default is 99 |
conf |
confidence level of the simulated polygons. Default is 0.95 |
diagfun |
user-defined function used to obtain the diagnostic measures from the fitted model object |
simfun |
user-defined function used to simulate a random sample from the model estimated parameters |
fitfun |
user-defined function used to re-fit the model to simulated data |
verb |
logical. If |
sort.res |
logical. If |
closest.angle |
logical. If |
angle.ref |
the reference angle from which points will be sorted starting from the closest angle to the input (in radians). Defaults to |
counter.clockwise |
logical. Should the points be ordered counter-clockwise or clockwise from the reference angle? |
xlab |
argument passed to |
ylab |
argument passed to |
main |
argument passed to |
clear.device |
logical. If |
point.col |
a vector of length 2 with the colors of the points that are inside and outside of the simulated polygons |
point.pch |
a vector of length 2 with the point characters of the points that are inside and outside of the simulated polygons |
... |
further arguments passed to |
x |
an object of class |
This approach relies on the same strategy used for producing half-normal plots with simulation envelopes. Given a vector of bivariate model diagnostics, the angle each point makes with the origin is calculated to order them. This can be fine-tuned using the logical arguments closest.angle
, angle.ref
, and counter.clockwise
, see the Arguments section above.
Then, sim
bivariate response variables are simulated from the fitted model, using the same model matrices, error distribution and fitted parameters, using the function defined as simfun
. The model is refitted to each simulated sample, using the function defined as fitfun
. Next, we obtain the same type of model diagnostics, using diagfun
, again ordered the same way the original bivariate sample was. We have, for each bivariate diagnostic, sim
simulated bivariate diagnostics forming the whole cloud of simulated diagnostics.
By default, we then obtain the convex hulls of each set of the $s$ sets of points and obtain a reduced polygon whose area is (conf * 100
)% of the original convex hull's area, forming the simulated polygon. This is equivalent to passing the argument reduce.polygon = "proportional"
to plot.bivrp
. The argument reduce.polygon = "bag"
can be used to obtain a (conf * 100
)% bagplot as the simulated polygon instead of a convex hull. The points are then connected to the centroids of their respective simulated polygons and, if they lie outside the polygons, they are drawn in red. For the final display, the polygons are erased so as to ease visualization.
There is no automatic implementation of a bivariate model in this function, and hence users must provide three functions for bivrp
. The first function, diagfun
, must extract the desired model diagnostics from a model fit object. The second function, simfun
, must return the response variable, simulated using the same error distributions and estimated parameters from the fitted model. The third and final function, fitfun
, must return a fitted model object. See the Examples section.
This function produces a plot by passing the computed object to plot.bivrp
. The print
method returns a data.frame
containing all ordered simulated bivariate diagnostics.
The function returns an object of class "bivrp", which is a list containing the following components:
reslist.ord |
list of ordered diagnostics from model refitting to each simulated dataset |
res.original.ord |
original model diagnostics |
res1 |
diagnostics from variable 1 |
res2 |
diagnostics from variable 2 |
res.original1 |
original model diagnostics for variable 1 |
res.original2 |
original model diagnostics for variable 2 |
conf |
confidence level of the simulated polygons |
Rafael A. Moral <rafael.deandrademoral@mu.ie> and John Hinde
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 | ## simulating a bivariate normal response variable
require(mvtnorm)
n <- 40
beta1 <- c(2, .4)
beta2 <- c(.2, .2)
x <- seq(1, 10, length = n)
X <- model.matrix(~ x)
mu1 <- X%*%beta1
mu2 <- X%*%beta2
sig1 <- 2
sig2 <- 3
sig12 <- -1.7
Sig1 <- diag(rep(sig1), n)
Sig2 <- diag(rep(sig2), n)
Sig12 <- diag(rep(sig12), n)
V <- rbind(cbind(Sig1, Sig12),
cbind(Sig12, Sig2))
set.seed(2016)
Y <- as.numeric(rmvnorm(1, c(mu1, mu2), V))
## code for fitting the model estimating covariance or not
bivnormfit <- function(Y, X, covariance) {
n <- nrow(X)
p <- ncol(X)
y <- cbind(Y[1:n],Y[(n+1):(2*n)])
XtXinv <- solve(crossprod(X, X))
beta.hat <- XtXinv %*% crossprod(X, y)
mu.hat <- X%*%beta.hat
sigma.hat <- 1/n * t(y - mu.hat) %*% (y - mu.hat)
if(!covariance) sigma.hat <- diag(diag(sigma.hat))
cov.betas <- sigma.hat %x% XtXinv
se.s1 <- sqrt(2*sigma.hat[1]^2/(n-p+1))
se.s2 <- sqrt(2*sigma.hat[4]^2/(n-p+1))
if(!covariance) se.s12 <- NA else {
rho <- sigma.hat[2]/sqrt(sigma.hat[1]*sigma.hat[4])
se.s12 <- sqrt((1+rho^2)*sigma.hat[1]*sigma.hat[4]/(n-p+1))
}
se.betas <- sqrt(diag(cov.betas))
se.sigma <- c(se.s1, se.s2, se.s12)
coefs <- c(beta.hat, sigma.hat[1], sigma.hat[4], sigma.hat[2])
names(coefs) <- c("beta1.0", "beta1.1", "beta2.0", "beta2.1", "sig1", "sig2", "sig12")
fitted <- c(mu.hat)
resid <- Y - fitted
Sig1 <- diag(rep(sigma.hat[1]), n)
Sig2 <- diag(rep(sigma.hat[4]), n)
Sig12 <- diag(rep(sigma.hat[2]), n)
V <- rbind(cbind(Sig1, Sig12),
cbind(Sig12, Sig2))
llik <- dmvnorm(Y, c(mu.hat), V, log = TRUE)
ret <- list("coefs" = coefs, "covariance" = covariance, "n" = n,
"X" = X, "fitted" = fitted, "resid" = resid, "loglik" = llik,
"Y" = Y, "se" = c(se.betas, se.sigma))
class(ret) <- "bivnormfit"
return(ret)
}
## fitting bivariate models with and without estimating covariance
fit0 <- bivnormfit(Y, X, covariance=FALSE)
fit1 <- bivnormfit(Y, X, covariance=TRUE)
## likelihood-ratio test
2*(fit0$loglik - fit1$loglik)
pchisq(54.24, 1, lower=FALSE)
## function for extracting diagnostics (raw residuals)
dfun <- function(obj) {
r <- obj$resid
n <- obj$n
return(list(r[1:n], r[(n+1):(2*n)]))
}
## function for simulating new response variables
sfun <- function(obj) {
n <- obj$n
fitted <- obj$fitted
sig1 <- obj$coefs[5]
sig2 <- obj$coefs[6]
if(obj$covariance) sig12 <- obj$coefs[7] else sig12 <- 0
Sig1 <- diag(rep(sig1), n)
Sig2 <- diag(rep(sig2), n)
Sig12 <- diag(rep(sig12), n)
V <- rbind(cbind(Sig1, Sig12),
cbind(Sig12, Sig2))
Y <- as.numeric(rmvnorm(1, c(mu1, mu2), V))
return(list(Y[1:n], Y[(n+1):(2*n)], "X" = obj$X,
"covariance" = obj$covariance))
}
## function for refitting the model to simulated data
ffun <- function(new.obj) {
Ynew <- c(new.obj[[1]], new.obj[[2]])
bivnormfit(Ynew, new.obj$X, new.obj$covariance)
}
## Bivariate residual plot for model 1 (without estimating covariance)
plot1 <- bivrp(fit0, diagfun=dfun, simfun=sfun, fitfun=ffun, verb=TRUE)
## without polygon area reduction
plot(plot1, conf=1)
## drawing polygons
plot(plot1, add.polygon=TRUE)
## without ordering
plot(plot1, theta.sort=FALSE, kernel=TRUE, add.dplots=TRUE, superpose=TRUE)
## Bivariate residual plot for model 2 (estimating covariance)
plot2 <- bivrp(fit1, diagfun=dfun, simfun=sfun, fitfun=ffun, verb=TRUE)
## without polygon area reduction
plot(plot2, conf=1)
## drawing polygons
plot(plot2, add.polygon=TRUE, conf=1)
## without ordering
plot(plot2, theta.sort=FALSE, kernel=TRUE, add.dplots=TRUE, superpose=TRUE)
|
Loading required package: MASS
Loading required package: mvtnorm
[1] -54.24361
[1] 1.774378e-13
Simulation 1 out of 99
Simulation 2 out of 99
Simulation 3 out of 99
Simulation 4 out of 99
Simulation 5 out of 99
Simulation 6 out of 99
Simulation 7 out of 99
Simulation 8 out of 99
Simulation 9 out of 99
Simulation 10 out of 99
Simulation 11 out of 99
Simulation 12 out of 99
Simulation 13 out of 99
Simulation 14 out of 99
Simulation 15 out of 99
Simulation 16 out of 99
Simulation 17 out of 99
Simulation 18 out of 99
Simulation 19 out of 99
Simulation 20 out of 99
Simulation 21 out of 99
Simulation 22 out of 99
Simulation 23 out of 99
Simulation 24 out of 99
Simulation 25 out of 99
Simulation 26 out of 99
Simulation 27 out of 99
Simulation 28 out of 99
Simulation 29 out of 99
Simulation 30 out of 99
Simulation 31 out of 99
Simulation 32 out of 99
Simulation 33 out of 99
Simulation 34 out of 99
Simulation 35 out of 99
Simulation 36 out of 99
Simulation 37 out of 99
Simulation 38 out of 99
Simulation 39 out of 99
Simulation 40 out of 99
Simulation 41 out of 99
Simulation 42 out of 99
Simulation 43 out of 99
Simulation 44 out of 99
Simulation 45 out of 99
Simulation 46 out of 99
Simulation 47 out of 99
Simulation 48 out of 99
Simulation 49 out of 99
Simulation 50 out of 99
Simulation 51 out of 99
Simulation 52 out of 99
Simulation 53 out of 99
Simulation 54 out of 99
Simulation 55 out of 99
Simulation 56 out of 99
Simulation 57 out of 99
Simulation 58 out of 99
Simulation 59 out of 99
Simulation 60 out of 99
Simulation 61 out of 99
Simulation 62 out of 99
Simulation 63 out of 99
Simulation 64 out of 99
Simulation 65 out of 99
Simulation 66 out of 99
Simulation 67 out of 99
Simulation 68 out of 99
Simulation 69 out of 99
Simulation 70 out of 99
Simulation 71 out of 99
Simulation 72 out of 99
Simulation 73 out of 99
Simulation 74 out of 99
Simulation 75 out of 99
Simulation 76 out of 99
Simulation 77 out of 99
Simulation 78 out of 99
Simulation 79 out of 99
Simulation 80 out of 99
Simulation 81 out of 99
Simulation 82 out of 99
Simulation 83 out of 99
Simulation 84 out of 99
Simulation 85 out of 99
Simulation 86 out of 99
Simulation 87 out of 99
Simulation 88 out of 99
Simulation 89 out of 99
Simulation 90 out of 99
Simulation 91 out of 99
Simulation 92 out of 99
Simulation 93 out of 99
Simulation 94 out of 99
Simulation 95 out of 99
Simulation 96 out of 99
Simulation 97 out of 99
Simulation 98 out of 99
Simulation 99 out of 99
18 out of 40 points out of polygons (45%).
13 out of 40 points out of polygons (32.5%).
18 out of 40 points out of polygons (45%).
Simulation 1 out of 99
Simulation 2 out of 99
Simulation 3 out of 99
Simulation 4 out of 99
Simulation 5 out of 99
Simulation 6 out of 99
Simulation 7 out of 99
Simulation 8 out of 99
Simulation 9 out of 99
Simulation 10 out of 99
Simulation 11 out of 99
Simulation 12 out of 99
Simulation 13 out of 99
Simulation 14 out of 99
Simulation 15 out of 99
Simulation 16 out of 99
Simulation 17 out of 99
Simulation 18 out of 99
Simulation 19 out of 99
Simulation 20 out of 99
Simulation 21 out of 99
Simulation 22 out of 99
Simulation 23 out of 99
Simulation 24 out of 99
Simulation 25 out of 99
Simulation 26 out of 99
Simulation 27 out of 99
Simulation 28 out of 99
Simulation 29 out of 99
Simulation 30 out of 99
Simulation 31 out of 99
Simulation 32 out of 99
Simulation 33 out of 99
Simulation 34 out of 99
Simulation 35 out of 99
Simulation 36 out of 99
Simulation 37 out of 99
Simulation 38 out of 99
Simulation 39 out of 99
Simulation 40 out of 99
Simulation 41 out of 99
Simulation 42 out of 99
Simulation 43 out of 99
Simulation 44 out of 99
Simulation 45 out of 99
Simulation 46 out of 99
Simulation 47 out of 99
Simulation 48 out of 99
Simulation 49 out of 99
Simulation 50 out of 99
Simulation 51 out of 99
Simulation 52 out of 99
Simulation 53 out of 99
Simulation 54 out of 99
Simulation 55 out of 99
Simulation 56 out of 99
Simulation 57 out of 99
Simulation 58 out of 99
Simulation 59 out of 99
Simulation 60 out of 99
Simulation 61 out of 99
Simulation 62 out of 99
Simulation 63 out of 99
Simulation 64 out of 99
Simulation 65 out of 99
Simulation 66 out of 99
Simulation 67 out of 99
Simulation 68 out of 99
Simulation 69 out of 99
Simulation 70 out of 99
Simulation 71 out of 99
Simulation 72 out of 99
Simulation 73 out of 99
Simulation 74 out of 99
Simulation 75 out of 99
Simulation 76 out of 99
Simulation 77 out of 99
Simulation 78 out of 99
Simulation 79 out of 99
Simulation 80 out of 99
Simulation 81 out of 99
Simulation 82 out of 99
Simulation 83 out of 99
Simulation 84 out of 99
Simulation 85 out of 99
Simulation 86 out of 99
Simulation 87 out of 99
Simulation 88 out of 99
Simulation 89 out of 99
Simulation 90 out of 99
Simulation 91 out of 99
Simulation 92 out of 99
Simulation 93 out of 99
Simulation 94 out of 99
Simulation 95 out of 99
Simulation 96 out of 99
Simulation 97 out of 99
Simulation 98 out of 99
Simulation 99 out of 99
10 out of 40 points out of polygons (25%).
8 out of 40 points out of polygons (20%).
8 out of 40 points out of polygons (20%).
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.