Nothing
#------------------------------------------------------#
# This function computes the Cov of std Reg coeffs for #
# fixed X using a method described by Yuan and Chan #
# (2011) #
# #
# #
# Arguments: #
# X - matrix of predictor scores #
# y - vector of criterion scores #
# cov.x - covariance matrix for predictors #
# cov.xy - vector of covariances between #
# predictors and criterion #
# var.y - criterion variance #
# var.error - var.error #
# Nobs - number of observations #
# #
# This function accepts either (1) raw data, or (2) #
# second-order moments (covariances or correlations) #
# and sample size. #
# #
# Output #
# seBeta - normal theory standard errors for #
# standarized regression coefficients #
# with Fixed predictors. #
# covBeta - normal theory covariance matrix of #
# standarized regsion coefficients #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#' Covariance Matrix and Standard Errors for Standardized Regression
#' Coefficients for Fixed Predictors
#'
#' Computes Normal Theory Covariance Matrix and Standard Errors for
#' Standardized Regression Coefficients for Fixed Predictors
#'
#'
#' @param X Matrix of predictor scores.
#' @param y Vector of criterion scores.
#' @param cov.x Covariance or correlation matrix of predictors.
#' @param cov.xy Vector of covariances or correlations between predictors and
#' criterion.
#' @param var.y Criterion variance.
#' @param var.error Optional argument to supply the error variance: var(y -
#' yhat).
#' @param Nobs Number of observations.
#' @return \item{cov.Beta}{Normal theory covariance matrix of standardized
#' regression coefficients for fixed predictors.} \item{se.Beta}{Standard
#' errors for standardized regression coefficients for fixed predictors.}
#' @author Jeff Jones and Niels Waller
#' @seealso \code{\link{seBeta}}
#' @references Yuan, K. & Chan, W. (2011). Biases and standard errors of
#' standardized regression coefficients. \emph{Psychometrika, 76(4)}, 670-690.
#' @keywords Statistics
#' @export
#' @import MASS
#' @examples
#'
#' ## We will generate some data and pretend that the Predictors are being held fixed
#'
#' library(MASS)
#' R <- matrix(.5, 3, 3); diag(R) <- 1
#' Beta <- c(.2, .3, .4)
#'
#' rm(list = ".Random.seed", envir = globalenv()); set.seed(123)
#' X <- mvrnorm(n = 200, mu = rep(0, 3), Sigma = R, empirical = TRUE)
#' y <- X %*% Beta + .64*scale(rnorm(200))
#'
#' seBetaFixed(X, y)
#'
#' # $covBeta
#' # b1 b2 b3
#' # b1 0.003275127 -0.001235665 -0.001274303
#' # b2 -0.001235665 0.003037100 -0.001491736
#' # b3 -0.001274303 -0.001491736 0.002830157
#' #
#' # $seBeta
#' # b1 b2 b3
#' # 0.05722872 0.05510989 0.05319922
#'
#' ## you can also supply covariances instead of raw data
#'
#' seBetaFixed(cov.x = cov(X), cov.xy = cov(X, y), var.y = var(y), Nobs = 200)
#'
#' # $covBeta
#' # b1 b2 b3
#' # b1 0.003275127 -0.001235665 -0.001274303
#' # b2 -0.001235665 0.003037100 -0.001491736
#' # b3 -0.001274303 -0.001491736 0.002830157
#' #
#' # $seBeta
#' # b1 b2 b3
#' # 0.05722872 0.05510989 0.05319922
#'
#'
seBetaFixed <- function(X=NULL,y=NULL,cov.x=NULL,cov.xy=NULL,
var.y=NULL,var.error = NULL,Nobs=NULL) {
########################################################################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~ Error Checking ~~~~~~~~~~~~~~~~~~~~~~~~~~~#
if(is.null(X) & !is.null(y))
stop("\n y is not defined\n Need to specify both X and y\n")
if(!is.null(X) & is.null(y))
stop("\n X is not defined\n Need to specify both X and y\n")
if(is.null(X) & is.null(y)) {
if(is.null(cov.x) | is.null(cov.xy) | is.null(var.y) | is.null(Nobs))
stop("\nYou need to specify all covariances and the number of observations\n")
N <- Nobs
p <- nrow(cov.x)
} else {
# if function is given data, compute all the covariances and variances
cov.x <- cov(X)
cov.xy <- cov(X,y)
var.y <- var(y)
N <- length(y)
p <- ncol(X)
}
# create diagonal matrix of predictor standard deviations (Dx),
# unstandardized regression coefficients (b), and
# error variance from regression model (var.error)
var.y <- as.numeric(var.y)
Dx <- diag(sqrt(diag(cov.x)))
b <- solve(cov.x)%*%cov.xy
# if error variance is not supplied, calculate it
if(is.null(var.error)) var.error <- as.numeric(var.y - t(b)%*%cov.x%*%b)*(N-1)/(N-p-1)
# create Gamma matrix which is the
# covariance matrix of b (unstandardized) and sigma^2 y (variance of y)
var.b <- solve(cov.x)*var.error
cov.b.var.y <- 2*b*var.error
var.var.y <- 4*(t(b)%*%cov.x%*%b)*var.error + 2*var.error^2
Gamma <- rbind(cbind(var.b,cov.b.var.y),c(cov.b.var.y,var.var.y))
# Create jacobian for the function g(b,var.y) = Dx%*%b/sd.y
deriv_b <- Dx*(1/sqrt(var.y))
deriv_var.y <- -Dx%*%b/(2*var.y^(3/2))
jacobian <- cbind(deriv_b, deriv_var.y)
# create covariance matrix of standardized regression coefficients
cov.mat <- jacobian%*%Gamma%*%t(jacobian)/N
# name the covariance matrix
b.nms <- NULL
for(i in 1:p) b.nms[i] <- paste("b",i,sep='')
rownames(cov.mat) <- colnames(cov.mat) <- b.nms
list(covBeta=cov.mat, seBeta=sqrt(diag(cov.mat)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.