Nothing
simulate.bic.glm <- function(object, nsim = 1,
seed = NULL, newdata=NULL,
type = c("coef", "link", "response"), ...){
##
## 1. seed?
##
# copy code from stats:::simulate.lm
if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1)
if (is.null(seed))
RNGstate <- get(".Random.seed", envir = .GlobalEnv)
else {
R.seed <- get(".Random.seed", envir = .GlobalEnv)
set.seed(seed)
RNGstate <- structure(seed, kind = as.list(RNGkind()))
on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
}
##
## 2. local copies of components of object
##
postprob <- object[['postprob']]
x <- object[['x']]
y <- object[['y']]
mle <- object[['mle']]
linkinv <- object[['linkinv']]
##
## 3. wt and family from call
##
cl <- as.list(object[['call']])
wt <- cl[['wt']]
fam <- cl[['glm.family']]
# convert fam to the form required by glm.fit
family <- get(as.character(fam),
mode = "function", envir = parent.frame())
if (is.function(family))
family <- family()
if (is.null(family$family)) {
print(family)
stop("'family' not recognized")
}
##
## 4. newdata?
##
if(is.null(newdata)){
newdata <- x
} else {
# ensure that newdata has the same variables as x
# Number of columns the same?
if(ncol(newdata) != ncol(x)){
colEr <- paste0("The columns of newdata ",
"must match object$x:\nncol(newdta) = ",
ncol(newdata), " != ncol(object$x) = ",
ncol(x) )
stop(colEr)
}
# Column names the same?
oopsNx <- which(names(newdata) != names(x))
if(length(oopsNx)>0){
# Does the names problem disappear
# by suppressing trailing ".x"?
namesx <- sub('\\.x$', '', names(x))
oopsNx2 <- which(names(newdata) != namesx)
if(length(oopsNx2)>0){
erMsg <- paste0(
"names(newdata) do not match those of",
" object$x.\nLook at column")
if(length(oopsNx)>1){
erMsg <- paste0(erMsg, "s")
}
stop(erMsg, " ", paste(oopsNx, collapse=", "))
}
}
# Factors in the same places?
Facs <- which(sapply(x, class)=="factor")
newFacs <- which(sapply(newdata, class)
== "factor")
if(length(Facs) != length(newFacs))
stop("newdata does not have the same ",
"number of factors as object$x: ",
"newdata has ", length(newFacs),
"; object$x has ", length(Facs))
# Factors have the same levels?
for(i in Facs){
lvlNi <- levels(newdata[[i]])
lvlxi <- levels(x[[i]])
if(length(lvlNi) != length(lvlxi))
stop("levels(newdata[[", i,
"]]) are different from levels(x[[",
i, "]])")
oopsLi <- which(lvlNi != lvlxi)
if(length(oopsLi)>0)
stop("levels(newdata[[", i,
"]]) are different from levels(x[[",
i, "]])")
}
}
##
## 5. model.matrices
##
xMat <- model.matrix(~., x)
newMat <- model.matrix(~., newdata)
##
## 6. sims = matrix to hold ouputs
##
nComponents <- length(postprob)
nobs <- NROW(newdata)
sims <- matrix(NA, nobs, nsim)
colnames(sims) <- paste0("sim_", 1:nsim)
if(!is.null(rownames(newMat)))
rownames(sims) <- rownames(newMat)
# if(is.null(wt))wt <- rep(1, nobs)
##
## 7. Which models where?
##
rmdl <- sample(1:nComponents, nsim, TRUE,
postprob)
#
nVars <- NCOL(mle)
Coefs <- matrix(0, nVars, nsim)
dimnames(Coefs) <- list(
colnames(newMat), colnames(sims) )
##
## 8. Simulate in link space
##
for(Comp in 1:nComponents){
nsimComp <- sum(rmdl==Comp)
if(any(rmdl==Comp)){
xMatC <- xMat[, mle[Comp,]!=0, drop=FALSE]
refitComp <- glm.fit(x=xMatC, y=y,
weights=wt, start=mle[Comp, mle[Comp,]!=0],
family=family)
class(refitComp) <- 'glm'
vc <- vcov(refitComp)
simCoef <- mvtnorm::rmvnorm(nsimComp,
coef(refitComp), vc)
Coefs[mle[Comp,]!=0, rmdl==Comp] <- t(simCoef)
#
newM <- newMat[, mle[Comp,]!=0, drop=FALSE]
predComp <- tcrossprod(newM, simCoef)
sims[, rmdl==Comp] <- predComp
}
}
##
## 9. return type(s) desired
##
# 9.1. only 1 type
if(length(type)<1)stop('No "type" requested.')
if(length(type)<2){
if(!is.na(pmatch(type, 'coef'))){
out <- Coefs
attr(out, 'seed') <- RNGstate
return(out)
} else if(!is.na(pmatch(type, 'link'))){
out <- data.frame(sims)
attr(out, 'seed') <- RNGstate
return(out)
} else if(is.na(pmatch(type, 'response'))){
stop('Not a recognized type. type = ', type)
} else {
out <- data.frame(linkinv(sims))
attr(out, 'seed') <- RNGstate
return(out)
}
}
# 9.2. more than 1 type requested
out <- vector('list', length(type))
names(out) <- type
if(any(!is.na(pmatch(type, 'coef')))){
out[['coef']] <- Coefs
}
if(any(!is.na(pmatch(type, 'link')))){
out[['link']] <- data.frame(sims)
}
if(any(!is.na(pmatch(type, 'response')))){
out[['response']] <- data.frame(linkinv(sims))
}
attr(out, 'seed') <- RNGstate
# cat('print(type)\n')
# print(type)
out
#
attr(out, 'seed') <- RNGstate
out
}
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.