# R/internal_functions.R In sValues: Measures of the Sturdiness of Regression Coefficients

#### Defines functions text_breakmelt.coefboundsbayesr2_combscombss_valuepriorVarmodel_frame_scale

```##### Internal Functions  - Not for user interaction #####

# scale data
##' @importFrom stats model.matrix
model_frame_scale <- function(formula, data, scale = TRUE){
x <- model.matrix(formula, data)[ ,-1, drop = FALSE]
if(scale) x  <- scale(x)
y <- data[deparse(formula[[2]])]
if(scale) y <- scale(y)
return(data.frame(y, x))
}

# Generates Prior Variance given favorites
.priorVar <- function(prior_R2, ols, favorites = NULL, R2_favorites = NULL){
k <- ols\$rank - 1
priors <- rep(prior_R2/(k), k)

if(!is.null(favorites)){
fav <- which(names(ols\$coefficients[-1] )%in% favorites)
if(length(fav) == 0) stop(paste("No variable matches", favorites))
if(length(fav) == length(names(ols\$coefficients[-1]))) stop("You can't choose all variables as favorites")
priors[fav]  <-  (R2_favorites)/(length(fav))
priors[-fav] <-  (prior_R2 - R2_favorites)/(length(priors[-fav]))
}

V <- diag(k)*(priors)

return(V)
}

# calculates s-values given a vector
.s_value <- function(vector){
min <- min(vector)
max <- max(vector)
(max +  min) / (max -  min)
}

# code originally from caTools, but caTools package is being retired by CRAN
combs <- function(v, k)
{
n = length(v)
if (n == k)
P = matrix(v, 1, n)
else if (k == 1)
P = matrix(v, n, 1)
else if (k == n - 1)
P = matrix(rep(v, each = n - 1), n, n - 1)
else if (k < n) {
P = matrix(0, 0, k)
if (k < n & k > 1) {
for (i in 1:(n - k + 1)) {
Q = combs(v[(i + 1):n], k - 1)
j = nrow(Q)
P = rbind(P, cbind(rep(v[i], j), Q))
}
}
}
else stop("combs: number m has to be smaller or equal to length of vector v")
return(P)
}

#
r2_combs <- function(R2_bounds){
r2_combs  <- as.data.frame(t(combs(R2_bounds, 2)))
ord       <- order(t(r2_combs)[,1, drop = FALSE], -order(t(r2_combs)[,2, drop = FALSE]))
r2_combs  <- r2_combs[,ord, drop = FALSE]
names(r2_combs) <- paste("R2", lapply(r2_combs, paste, collapse = "_"), sep = "_")
r2_combs
}

# Calculates the bayesian estimates given a prior R2
##' @importFrom stats vcov
.bayes <- function(ols, prior_R2, favorites = NULL, R2_favorites = NULL){

if(class(ols)!= "lm") stop("ols must be an object of class 'lm'")

#--- prior variance
V <- .priorVar(prior_R2, ols, favorites, R2_favorites)
#--- ols coefficients

b <- ols\$coefficients[-1]

#--- ols precision
H <- solve(vcov(ols)); H <- H[-1, , drop = FALSE]; H <- H[,-1, drop = FALSE]

#-- bayes estimate
bayes <- solve((H+solve(V)))%*%H%*%b

colnames(bayes) <- paste("b_bayes", prior_R2, sep="_")

#-- bayes var-covar
bayes_cov <- solve((H+solve(V)))

#-- bayes t
t <- bayes/sqrt(diag(bayes_cov))
colnames(t) <- paste("t_bayes", prior_R2, sep="_")

#-- formatting data before returning
coefs <- data.frame(coef=bayes)
names(coefs) <- paste(c("b"), prior_R2, sep="_")

ret <- list(coefficients=bayes, vcov=bayes_cov,  t=t)
return(ret)
}

# Extreme bounds given R2 bounds
##' @importFrom stats vcov
.bounds <- function(ols, coef_indices, R2_bounds, favorites = NULL, R2_favorites = NULL){

#-- linear combination --#
k <- ols\$rank - 1
phi <- rep(0, k)
phi[coef_indices] <- 1

#-- prior variances --#
V1 <- .priorVar(min(R2_bounds), ols, favorites, min(R2_favorites))
V2 <- .priorVar(max(R2_bounds), ols, favorites, max(R2_favorites))

#--- Extreme bounds --#
b <- ols\$coefficients[-1]
H <- solve(vcov(ols)); H <- H[-1, , drop = FALSE]; H <- H[,-1, drop = FALSE]
G <- (H + solve(V2))%*%solve((solve(V1)-solve(V2)))%*%(H+solve(V2)) + (H+solve(V2))
f <- solve(H+solve(V1))%*%(H%*%b + (solve(V1)-solve(V2))%*%solve(H+solve(V2))%*%(H%*%b)/2)
c <- t(b)%*%H%*%solve(H + solve(V2))%*%(solve(V1)-solve(V2))%*%solve(H+solve(V1))%*%(H%*%b)/4
up <- t(phi)%*%f + sqrt(t(phi)%*%solve(G)%*%phi)%*%sqrt(c)
low <- t(phi)%*%f - sqrt(t(phi)%*%solve(G)%*%phi)%*%sqrt(c)
data.frame(low, up)
}

# helper function to melt coefficients for plotting
.melt.coef <- function(x, type){
coef <- coef(x, type)
coef\$coefficients <- rownames(coef)
melted <- suppressMessages(melt(coef, stringsAsFactors=FALSE, value.name = type))
return(melted)
}

# breaks text for pretty printing
text_break <- function(text, print.length = 5){
if(print.length < 1) print.length <- 1
if(length(text) > print.length){
more <- paste("and", length(text) - print.length, "more.")
text <- c(text[1:print.length], more)
}
text
}
```

## Try the sValues package in your browser

Any scripts or data that you put into this service are public.

sValues documentation built on July 16, 2018, 1:04 a.m.