# R/predictWorkers.R In texmex: Statistical Modelling of Extreme Values

```addCov <- function(res, X){ # used in linearPredictors.* to add covariates to columns reported in output
if(!is.null(dim(X))){
if(dim(X)[2] > 1){
cov <- X[,colnames(X) != "(Intercept)", drop=FALSE]
res <- cbind(res, cov)
if(is.vector(cov)) colnames(res)[dim(res)[2]] <- colnames(X)[colnames(X) != "(Intercept)"]
}
else {
if( any(X != 1) ){
res <- cbind(res, X)
}
}
}
res
}

texmexMakeParams <-
# Take parameter vector and list of datasets to compute
# parameters for each row of the data.
function(co, data){

np <- length(data)
p <- vector('list', length=np)

wh <- 1
for (i in 1:np){
which <- wh:(wh -1 + ncol(data[[i]]))
p[[i]] <- c(co[which] %*% t(data[[i]]))
wh <- wh + ncol(data[[i]])
}

do.call('cbind', p)
}

texmexMakeCovariance <-
# Get covariance matrix for each parameter and get coefficents for each row of the data
function(cov, data){
covs <- v <- vector('list', length=length(data))

# First get covariance matrix for each main parameter (e.g. phi=a'X1, xi=b'X2)
wh <- 1
for (i in 1:length(covs)){
which <- wh:(wh -1 + ncol(data[[i]]))
covs[[i]] <- cov[which, which, drop=FALSE]
# covs[[]i] contains the block of the full covariance which relates to parameter[i]

# Get the variance of the linear predictors
v[[i]] <- rowSums((data[[i]] %*% covs[[i]]) * data[[i]])

wh <- wh + ncol(data[[i]])
}
names(v) <- names(D)
# Each element of v contains the variance for one linear predictor for every observation

# We now need the off-diagonal elements of the covariance matrix. The dimensions
# of the covariance will depend on the length of object\$data\$D

getOffDiagonal <- function(row, x1, x2, cov, cindex){
covar <- 0
for (j in 1:ncol(x1)){
for (k in 1:ncol(x2)){
covar <- covar + x1[row, j] * x2[row, k] * cov[j, cindex + k] # XXX <-----------------------------------
} # Close for j
} # Close for i
covar
} # Close getOffDiagonal

getCovEntry <- function(cov, data){
# Recursively produce off-diagonal elements of the covariance.
# These are computed by row.
x1 <- data[[1]]
data[[1]] <- NULL
n <- length(data)

res <- vector('list', length=n)
wh <- 1
for (i in 1:n){
res[[i]] <- sapply(1:nrow(x1), getOffDiagonal, x1, data[[i]], cov, wh)
wh <- wh + ncol(data[[i]]) # Increment column index for next dataset
}
if (n > 1){
# Get next covariance block. Already pruned data, now prune cov of those elements already used.
cov <- cov[-(1:ncol(x1)), -(1:ncol(x1)), drop=FALSE]
res <- c(res, getCovEntry(cov, data))
}
else {
res
}
} # Close getCovEntry

co <- getCovEntry(cov, data)

# Now need to restructure to return a list of covariance matrices,
# one element for every observation. The covariance is of the linear
# predictors.
getM <- function(i, variance, covariance){
va <- lapply(variance, function(x){ x[[i]] })
co <- lapply(covariance, function(x){ x[[i]] })

res <- diag(unlist(va))
res[upper.tri(res)] <-  res[lower.tri(res)] <- unlist(co)
res
}

lapply(1:nrow(data[[1]]), getM, v, co)
}

texmexPredictSE <- function(x){
# x is an object returned by texmexMakeCovariance. It
# is a list, each element of which is a covariance matrix.

se <- sapply(x, function(x){ sqrt(diag(x)) })
}

texmexMakeCI <-
# Compute CIs from point estimates and standard errors
function(params, ses, alpha){
z <- qnorm(1 - alpha/2)

ci <- lapply(1:ncol(params), function(i, x, s, z){
cbind(x[, i] - z*s[, i], x[, i] + z*s[, i])
}, x=params, s=ses, z=z)
ci <- do.call('cbind', ci)

nms <- rep('', ncol(ci))
nms[(1:ncol(ci)) %% 2 == 1] <- paste(colnames(params), '.lo', sep = '')
nms[(1:ncol(ci)) %% 2 == 0] <- paste(colnames(params), '.hi', sep = '')
colnames(ci) <- nms

ci
}

if (is.null(newdata)){
res <- x\$data\$D
} else {
xl <- function(i, fo, data, xlev){
if (length(xlev[[i]]) > 0){
model.matrix(fo[[i]], data, xlev=xlev[[i]][[1]])
} else {
model.matrix(fo[[i]], data)
}
} # Close function
fo <- x\$formulae
res <- lapply(1:length(fo), xl, fo=fo, data=newdata, xlev=x\$xlevels)
names(res) <- names(fo)
}
invisible(res)
}

texmexMakeCISim <- function(x, alpha, object, sumfun, M){
if (is.null(sumfun)){
sumfun <- function(x){
c(Mean=mean(x), quantile(x, prob=c(.50, alpha/2,  1 - alpha/2)) )
}
neednames <- TRUE
} else {
neednames <- FALSE
}

t(apply(x, 2, sumfun))
}
```

## Try the texmex package in your browser

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

texmex documentation built on May 2, 2019, 5:41 a.m.