Nothing
## SimDat models
## used for simulation of random variables
setClass("Model")
#setClass("GlmModel",
# contains="Model",
# representation(
# formula="formula",
# family="family",
# coefficients="numeric",
# parameters=list() # nuisance parameters (e.g., sd)
# )
#)
setClass("ParModel",
representation(
formula="formula",
coefficients="vector" # vector or (column) matrix
)
)
setMethod("isRandom","ParModel",
function(object,...) {
FALSE
}
)
# TODO: use formula with conditional instead of grouping
setClass("RandomParModel",
contains="ParModel",
representation(
#random.formula="formula", # formula specifying random effects models
#random.coefficients = "matrix", # matrix with random coefficients
grouping="numeric", # index variable for grouping levels (repetitions) of random coefficients
sigma="array" # ncol(coeff) * ncol(coeff) * (1 | unique(grouping)) array with covariance matrices for random coefficients
)
)
setMethod("isRandom","RandomParModel",
function(object,...) {
TRUE
}
)
#setMethod("simulate",signature(object="Variable",model="RandomParModel"),
# function(object,nsim=1,seed,model,...) {
# object
# }
#)
setMethod("simulate",signature(object="RandomParModel"),
function(object,nsim=1,seed,...) {
ns <- dim(object@sigma)[3]
nc <- dim(object@sigma)[1]
if(!is.matrix(object@coefficients)) {
if(length(object@coefficients) == nc) {
ng <- 1
if(nc == 1) {
coeff <- rnorm(ng,mean=0,sd=as.numeric(sqrt(object@sigma[,,1])))
} else {
coeff <- as.numeric(rmvnorm(ng,mean=rep(0,nc),sigma=object@sigma[,,1]))
}
} else {
if(nc == 1) {
ng <- length(object@coefficients)
coeff <- rnorm(ng,mean=0,sd=as.numeric(sqrt(object@sigma[1,1,])))
} else {
stop("coefficients doesn't seem to have the correct size")
}
}
} else {
coeff <- object@coefficients
ng <- nrow(coeff)
if(nc > 1) {
if(ns > 1) {
if(ns == ng) {
for(i in 1:ng) {
coeff[i,] <- rmvnorm(1,mean=rep(0,nc),sigma=object@sigma[,,i])
}
} else {
stop("size of sigma doesn't match grouping")
}
} else {
coeff <- rmvnorm(ng,mean=rep(0,nc),sigma=object@sigma[,,1])
}
} else {
if(ncol(coeff) == 1) {
coeff <- matrix(rnorm(ng,mean=0,sd=sqrt(object@sigma[1,1,])),ncol=1)
} else {
stop("sigma is a 1 by 1 matrix, but coefficients has more than 1 column")
}
}
}
object@coefficients <- coeff
object
}
)
setClass("MixedParModel",
contains="ParModel",
representation(
random="RandomParModel"
)
)
setMethod("simulate",signature(object="MixedParModel"),
function(object,nsim=1,seed,...) {
call <- match.call()
args <- as.list(call[2:length(call)])
args[["object"]] <- object@random
object@random <- do.call("simulate",args=args)
object
}
)
setMethod("isRandom","MixedParModel",
function(object,...) {
TRUE
}
)
setClass("RandomParModelList",
contains="list"
)
setMethod("isRandom","RandomParModelList",
function(object,...) {
any(unlist(lapply(object,isRandom)))
}
)
RandomParModelList <- function(list) {
if(!is.list(list)) stop("argument is not a list")
else {
out <- new("RandomParModelList",
.Data <- list)
}
out
}
setValidity("RandomParModelList",
method=function(object) {
any(unlist(lapply(object,function(x) is(x,"RandomParModel"))==TRUE))
}
)
setClass("MultilevelParModel",
contains="ParModel",
representation(
random="RandomParModelList"
)
)
setMethod("isRandom","MultilevelParModel",
function(object,...) {
isRandom(object@random,...)
}
)
setClass("VecParModel",
contains="ParModel"
)
setClass("MatParModel",
contains="ParModel"
)
setMethod("predict",signature(object="ParModel"),
function(object,...,data,na.action = c("exclude","pass")) {
na.action <- match.arg(na.action)
if(!missing(data)) {
m <- model.matrix(object@formula,...,data=data)
} else {
if(length(attr(terms(object@formula),"order")) == 0 & attr(terms(object@formula),"response") == 0) {
# intercept only model
m <- matrix(1,ncol=1,nrow=1)
} else {
m <- model.matrix(object@formula,...)
}
}
#TODO: don't know if this is general enough!
coeff <- object@coefficients
excl <- is.na(coeff)
if(sum(excl) > 0 & na.action=="exclude") {
m <- m[,!excl]
coeff <- coeff[!excl]
}
as.numeric(m%*%coeff)
}
)
setMethod("predict",signature(object="RandomParModel"),
function(object,...,data) {
if(!missing(data)) {
m <- model.matrix(object@formula,data=data)
} else {
m <- model.matrix(object@formula)
}
if(is.matrix(object@coefficients)) out <- rowSums(m*object@coefficients[object@grouping,]) else out <- rowSums(m*object@coefficients[object@grouping])
out
}
)
setMethod("predict",signature(object="MixedParModel"),
function(object,...,data) {
#call <- match.call()
#args[[1]] <- NULL
#args <- as.list(args)
#args <- as.list(call[2:length(call)])
#fp <- do.call("callNextMethod",args=args)
if(!missing(data)) {
fp <- callNextMethod(object=object,data=data,...)
rp <- predict(object@random,data=data,...)
} else {
fp <- callNextMethod(object=object,...)
rp <- predict(object@random,...)
}
#rp
rp + fp
#call
}
)
setMethod("predict",signature(object="VecParModel"),
function(object,...,data) {
if(!missing(data)) {
m <- model.matrix(object@formula,data=data)
} else {
m <- model.matrix(object@formula)
}
m%*%(object@coefficients)
}
)
setOldClass("gamlss.family")
setClass("GamlssModel",
contains="Model",
representation(
mu="ParModel",
sigma="ParModel",
nu="ParModel",
tau="ParModel",
family="gamlss.family"
)
)
setMethod("summary",signature(object="GamlssModel"),
function(object,data,type=c("aov","lm","glm"),SStype=3,...) {
type <- match.arg(type)
if(length(object@mu@formula) > 2) {
if(type=="aov") {
if(SStype == 1) {
return(summary(aov(formula=object@mu@formula,data=data,...)))
} else {
mod <- lm(formula=object@mu@formula,data=data,...)
return(Anova(mod,type=SStype,singular.ok=TRUE,...))
}
}
if(type=="glm") {
return(summary(glm(formula=object@mu@formula,data=data,...)))
}
if(type=="lm") {
return(summary(lm(formula=object@mu@formula,data=data,...)))
}
} else {
return(NULL)
}
}
)
setMethod("simulateFromModel",signature(object="RandomVariable",model="MixedParModel"),
function(object,model,nsim=1,seed,...,data) {
call <- match.call()
args <- as.list(call[2:length(call)])
DV <- do.call("simulate",args=args)
DV
}
)
# TODO: for consistency, should really simulate a DV, with a model as argument
setMethod("simulateFromModel",signature(object="RandomVariable",model="GamlssModel"),
function(object,model,nsim=1,seed,...,data) {
#if(isRandom(model))
model <- simulate(model,...)
npar <- model@family$nopar
args <- list()
if(!missing(data)) {
if(npar > 3) args[["tau"]] <- model@family$tau.linkinv(predict(model@tau,data=data))
if(npar > 2) args[["nu"]] <- model@family$nu.linkinv(predict(model@nu,data=data))
if(npar > 1) args[["sigma"]] <- model@family$sigma.linkinv(predict(model@sigma,data=data))
args[["mu"]] <- model@family$mu.linkinv(predict(model@mu,data=data))
} else {
if(npar > 3) args[["tau"]] <- model@family$tau.linkinv(predict(model@tau))
if(npar > 2) args[["nu"]] <- model@family$nu.linkinv(predict(model@nu))
if(npar > 1) args[["sigma"]] <- model@family$sigma.linkinv(predict(model@sigma))
args[["mu"]] <- model@family$mu.linkinv(predict(model@mu))
}
args[["n"]] <- length(object) #else args[["n"]] <- length(args$mu)
#if(!missing(DV)) {
if(isMetric(object)) {
min <- min(object)
max <- max(object)
w <- which(abs(c(min,max)) < Inf)
if(length(w) > 0) {
# have to truncate
if(length(w) == 2) {
type <- "both"
range <- c(min,max)
} else {
if(w == 1) {
type <- "left"
range <- min
} else {
type <- "right"
range <- max
}
}
rfun <- trun.r(range,family=model@family$family[1],type=type)
out <- do.call(rfun,args=args)
} else {
# no truncation
out <- do.call(paste("r",model@family$family[1],sep=""),args=args)
}
if(length(object@digits) > 0) out <- round(out,digits=object@digits)
} else {
# no truncation for non-metric variables
out <- do.call(paste("r",model@family$family[1],sep=""),args=args)
}
object@.Data <- out
#} else {
# DV <- do.call(paste("r",object@family$family[1],sep=""),args=args)
#}
return(object)
}
)
setClass("GammlssModel",
contains="GamlssModel"#,
)
setMethod("simulate",signature(object="GamlssModel"),
function(object,nsim=1,seed,...) {
return(object)
}
)
setMethod("simulate",signature(object="GammlssModel"),
function(object,nsim=1,seed,...) {
# simulate random coefficients
call <- match.call()
args <- as.list(call[2:length(call)])
args[["object"]] <- NULL
#targs <- args
#targs[["object"]] <- NULL
if(isRandom(object@mu)) object@mu <- do.call("simulate",args=c(object=object@mu,args))
if(isRandom(object@sigma)) object@sigma <- do.call("simulate",args=c(object=object@sigma,args))
if(isRandom(object@nu)) object@nu <- do.call("simulate",args=c(object=object@nu,args))
if(isRandom(object@tau)) object@tau <- do.call("simulate",args=c(object=object@tau,args))
object
#DV <- do.call("callNextMethod",args=args)
#DV
}
)
setClass("MvnormModel",
contains="Model",
representation(
mu="VecParModel",
sigma="array",
sigmaID="numeric"
)
)
setMethod("simulateFromModel",signature(object="RandomVariable",model="MvnormModel"),
function(object,model,nsim=1,seed,...,DV,data) {
if(!missing(data)) mu <- predict(object@mu,data=data,...) else mu <- predict(object@mu,...)
if(is(DV,"VariableList")) {
nDV <- length(DV)
# check for truncation
#trunc <- FALSE
trunc.points <- matrix(NA,ncol=2,nrow=length(DV))
for(i in 1:nDV) {
trunc.points[i,] <- c(min(DV[[i]]),max(DV[[i]]))
}
if(any(abs(trunc.points) < Inf)) {
# need a truncated MV normal
for(i in unique(object@sigmaID)) out[object@sigmaID==i,] <- rtmvnorm(nrow(mu),mean=rep(0,ncol(mean)),sigma=object@sigma[,,i])
} else {
for(i in unique(object@sigmaID)) out[object@sigmaID==i,] <- rmvnorm(nrow(mu),mean=rep(0,ncol(mean)),sigma=object@sigma[,,i])
}
} else {
if(is(DV,"Variable")) {
} else {
stop("DV should be a VariableList or Variable")
}
}
}
)
setClass("UniformModel",
contains = "GamlssModel"
)
UniformModel <- function() {
fam <- UN()
new("UniformModel",
family=fam)
}
setMethod("simulateFromModel",signature(object="RandomVariable",model="UniformModel"),
function(object,model,nsim=1,seed,...) {
n <- length(object)
if(isMetric(object)) {
min <- min(object)
max <- max(object)
if(min >= max) stop("maximum should be larger than minimum")
if(abs(min) == Inf) min <- -1e308
if(abs(max) == Inf) max <- 1e308
out <- runif(n,min,max)
if(length(object@digits) > 0) out <- round(out,digits=object@digits)
object@.Data <- out
} else {
levs <- levels(as.factor(object))
object@.Data <- sample(unique(object@.Data),size=n,replace=TRUE)
}
return(object)
}
)
#setClass("NormalModel",
# contains = "Model",
# representation(
# mean="numeric",
# sd="numeric"
# )
#)
NormalModel <- function(mean=0.0,sd=1.0) {
fam <- NO()
new("GamlssModel",
mu=new("ParModel",
formula=~1,
coefficients=fam$mu.linkfun(mean)),
sigma=new("ParModel",
formula=~1,
coefficients=fam$sigma.linkfun(sd)),
family=fam)
}
setClass("GamlssModel",
contains="Model",
representation(
mu="ParModel",
sigma="ParModel",
nu="ParModel",
tau="ParModel",
family="gamlss.family"
)
)
#setMethod("simulateFromModel",signature(object="RandomVariable",model="NormalModel"),
# function(object,model,nsim=1,seed,...) {
# n <- length(object)
# family <- NO()
# args <- list(n=n,mu=model@mean,sigma=model@sd)
# if(isMetric(object)) {
# min <- min(object)
# max <- max(object)
# w <- which(abs(c(min,max)) < Inf)
# if(length(w) > 0) {
# # have to truncate
# if(length(w) == 2) {
# type <- "both"
# range <- c(min,max)
# } else {
# if(w == 1) {
# type <- "left"
# range <- min
# } else {
# type <- "right"
# range <- max
# }
# }
# rfun <- trun.r(range,family=family$family[1],type=type)
# out <- do.call(rfun,args=args)
# } else {
# # no truncation
# out <- do.call(paste("r",family$family[1],sep=""),args=args)
# }
# #out <- rnorm(n,mean=model@mean,sd=model@sd)
# if(length(object@digits) > 0) out <- round(out,digits=object@digits)
# object@.Data <- out
# } else {
# stop("Cannot simulate a non-metric Variable from a NormalModel")
# }
# return(object)
# }
#)
setClass("ModelList",
contains="list"
)
ModelList <- function(list) {
if(!is.list(list)) stop("argument is not a list")
else {
out <- new("ModelList",
.Data <- list)
}
out
}
setValidity("ModelList",
method=function(object) {
all(lapply(object,function(x) is(x,"Model"))==TRUE)
}
)
setMethod("mean",signature(x="GamlssModel"),
function(x,data,...) {
if(!missing(data)) {
mu <- x@family$mu.linkinv(predict(x@mu,data=data))
} else {
mu <- x@family$mu.linkinv(predict(x@mu))
}
return(mu)
}
)
setMethod("sd",signature(x="GamlssModel"),
function(x,data,...) {
if(!missing(data)) {
sigma <- x@family$sigma.linkinv(predict(x@sigma,data=data))
} else {
sigma <- x@family$sigma.linkinv(predict(x@sigma))
}
return(sigma)
}
)
#setMethod("mean",signature(x="NormalModel"),
# function(x,...) {
# return(x@mean)
# }
#)
setMethod("mean",signature(x="UniformModel"),
function(x,DV,...) {
if(!is(DV,"RandomVariable") || !isMetric(DV)) stop("cannot compute the mean of a non-random and/or non metric Variable")
return(mean(c(min(DV),max(DV))))
}
)
setMethod("sd",signature(x="UniformModel"),
function(x,DV,...) {
if(!is(DV,"RandomVariable") || !isMetric(DV)) stop("cannot compute the mean of a non-random and/or non metric Variable")
return(sqrt((1/12)*(max(DV)-min(DV))^2))
}
)
#setMethod("sd",signature(x="NormalModel"),
# function(x,...) {
# return(x@sd)
# }
#)
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.