Nothing
# Copyright (C) 2011-2014 Thomas W. D. Möbius (kontakt@thomasmoebius.de)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see
# <http://www.gnu.org/licenses/>.
#' Design: Gaussian responses (known heteroscedasticity)
#'
#' Method for generating a sampling design for data generation
#' following a random effects meta regression model with
#' known heteroscedasticity.
#'
#' Generates a sampling design for the heterogeneity 'h' and
#' a heteroscedasticity 'd1', ..., 'dk'.
#'
#' Points in the design are selected via a maxi-min hypercube sampling
#' using the 'lhs' package in a predefined parameter cube.
#'
#' @param n resolution of the heterogeneity and heteroscedasticity
#' parameters, i.e. the number of of different (heterogeneity,
#' heteroscedasticity) pairs in the design.
#' @param h_bounds bounds of the heterogeneity.
#' @param d_bounds bounds of the heteroscedasticity.
#' @param x design matrix.
#' @return
#' Function returns a data frame. Each line of this data frame
#' can be an input to the function 'rY' which is used to sample
#' data from such a design.
#' @examples
#' dY <- designY(n=15L, h_bounds=c(0,1), d_bounds=c(0.01,2),
#' x=cbind(1,1:7))
#'
#' if(!all(dim(dY) == c(15,dim(cbind(1,1:7))[1]+1))) {
#' stop("Wrong dimension")
#' }
#' @export
designY <- function( n # resolution
, h_bounds # bounds of the heterogeneity
, d_bounds # bounds of the heteroscedasticity
, x # design matrix
) {
checkArg(n, "integer", len=1, lower=1, na.ok=FALSE)
checkArg(h_bounds, "numeric", len=2, lower=0, na.ok=FALSE)
checkArg(d_bounds, "numeric", len=2, lower=0, na.ok=FALSE)
ps <- makeParamSet( makeNumericParam("h" , lower=h_bounds[1],
upper=h_bounds[2])
, makeNumericVectorParam( "d"
, len=dim(x)[1]
, lower=d_bounds[1]
, upper=d_bounds[2]))
return(generateDesign(n=n, ps, fun=maximinLHS))
}
#' Design: Gaussian responses (unknown heteroscedasticity)
#'
#' Method for generating a sampling design for data generation
#' following a random effects meta regression model with
#' unknown heteroscedasticity.
#'
#' Generates a sampling design for the heterogeneity 'h',
#' heteroscedasticity 'd1', ..., 'dk', and study sizes 's1', ..., 'sk'.
#' This design can be used for testing methods that adjust for
#' uncertainty in the heteroscedasticity estimates by additionally
#' considering the size of the respected studies.
#'
#' Points in the design are selected via a maxi-min hypercube sampling
#' using the 'lhs' package in a predefined parameter cube.
#'
#' @param n resolution of the heterogeneity and heteroscedasticity
#' parameters, i.e., the number of of different (heterogeneity,
#' heteroscedasticity, sizes) tuple in the design.
#' @param h_bounds bounds of the heterogeneity.
#' @param d_bounds bounds of the heteroscedasticity.
#' @param s_bounds bounds of the study sizes.
#' @param x design matrix.
#' @return
#' Function returns a data frame. Each line of this data frame
#' can be an input to the function 'rD' which is used to sample
#' data from such a design.
#' @examples
#' dD <- designD(n=15L, h_bounds=c(0,1), d_bounds=c(0.01,2),
#' s_bounds=c(200L,2000L), x=cbind(1,1:7))
#'
#' if(!all(dim(dD) == c(15,2*dim(cbind(1,1:7))[1]+1))) {
#' stop("Wrong dimension")
#' }
#' @export
designD <- function( n # resolution
, h_bounds # bounds of the heterogeneity
, d_bounds # bounds of the heteroscedasticity
, s_bounds # bounds of the study sizes
, x # design matrix
) {
checkArg(n, "integer", len=1, lower=1, na.ok=FALSE)
checkArg(h_bounds, "numeric", len=2, lower=0, na.ok=FALSE)
checkArg(d_bounds, "numeric", len=2, lower=0, na.ok=FALSE)
checkArg(s_bounds, "integer", len=2, lower=1, na.ok=FALSE)
ps <- makeParamSet( makeNumericParam("h" , lower=h_bounds[1],
upper=h_bounds[2])
, makeNumericVectorParam( "d"
, len=dim(x)[1]
, lower=d_bounds[1]
, upper=d_bounds[2]))
psSizes <- makeParamSet(makeIntegerVectorParam( "s"
, len=dim(x)[1]
, lower=s_bounds[1]
, upper=s_bounds[2]))
return(cbind( generateDesign(n=n, ps, fun=maximinLHS)
, generateDesign(n=n, psSizes, fun=maximinLHS)))
}
#' Design: Binomial responses
#'
#' Method for generating a sampling design for data generation
#' following a binomial-Gaussian model.
#'
#' Generates a sampling design for the heterogeneity 'h', balancing
#' factors 'a1', ..., 'ak' of group assignments, and study sizes 's1',
#' ..., 'sk'. This design can be used for testing methods for inference
#' for the random effects meta regression model since the logarithm of
#' relative risks of each study is approximately Gaussian distributed.
#' One may use methods that adjust for uncertainty in the
#' heteroscedasticity estimates by additionally considering the size of
#' the respected studies.
#'
#' Points in the design are selected via a maxi-min hypercube sampling
#' using the 'lhs' package in a predefined parameter cube.
#'
#' @param n resolution of the heterogeneity. n is the number of
#' of different heterogeneity parameters in the design.
#' @param h_bounds bounds of the heterogeneity.
#' @param a_bounds bounds of the balancing factor of group assignments.
#' @param s_bounds bounds of the study sizes.
#' @param r fixed risk in the control.
#' @param x design matrix.
#' @return
#' Function returns a data frame. Each line of this data frame
#' can be an input to the function 'rB' which is used to sample
#' data from such a design.
#' @examples
#' dB <- designB(n=15L, h_bounds=c(0,1), a_bounds=c(-.3,3),
#' s_bounds=c(200L,2000L), r=0.03, x=cbind(1,1:5))
#'
#' if(!all(dim(dB) == c(15,2*dim(cbind(1,1:5))[1]+2))) {
#' stop("Wrong dimension")
#' }
#' @export
designB <- function( n
, h_bounds
, a_bounds
, s_bounds
, r
, x
) {
checkArg(n, "integer", len=1, lower=1, na.ok=FALSE)
checkArg(h_bounds, "numeric", len=2, lower=0, na.ok=FALSE)
checkArg(a_bounds, "numeric", len=2, na.ok=FALSE)
checkArg(s_bounds, "integer", len=2, lower=1, na.ok=FALSE)
checkArg(r, "numeric", len=1, lower=.0001, upper=.9999)
ps <- makeParamSet( makeNumericParam("h"
, lower=h_bounds[1]
, upper=h_bounds[2])
, makeNumericVectorParam( "a"
, len=dim(x)[1]
, lower=a_bounds[1]
, upper=a_bounds[2]))
psSizes <- makeParamSet(makeIntegerVectorParam( "s"
, len=dim(x)[1]
, lower=s_bounds[1]
, upper=s_bounds[2]))
return(cbind( generateDesign(n=n, ps, fun=maximinLHS)
, generateDesign(n=n, psSizes, fun=maximinLHS)
, data.frame(r=r)))
}
#' Data generation: Sampling data of clinical trials
#'
#' A random draw of a hierarchical binomial Gaussian model.
#'
#' It is always assumed that at least one response in a study has
#' happend, i.e., a response of 0 in a treatment or control group is
#' rounded up to 1. Note that this may lead to an overestimation of
#' small risks. If possible, make sure your sample sizes are large
#' enough to compensate for this effect.
#'
#' You may work around this by increasing study sizes.
#' @param h heterogeneity.
#' @param s study sizes.
#' @param a balance of group assignments.
#' @param r fixed risk in the control group.
#' @param x design matrix.
#' @param b regression coefficients.
#' @return A list containing the risk and a data frame with the studies.
#' @examples
#' h_test <- .03
#' s_test <- rep(2000, 13)
#' a_test <- rep(.3, 13)
#' x_test <- cbind(1,1:13)
#' b_test <- c(0.02, 0.03)
#' dat <- rBinomGauss(h=h_test, s=s_test, a=a_test, r=0.03 , x=x_test,
#' b=b_test)$study
#'
#' if(!all(dim(dat) == c(dim(x_test)[1], 4))) stop("Wrong dimension")
#' @export
rBinomGauss <- function( h # heterogeneity
, s # study sizes
, a # balance of group assignments
, r # fixed risk in the control group
, x # design matrix
, b # regression coefficients
) {
k <- length(s)
y <- rnorm(k, mean=as.vector(x%*%b), sd=sqrt(h))
p <- max(min(r * exp(-y), 0.999), 0.001)
sizeP <- round(s*(a+1)/2)
sizeR <- s - sizeP
n1 <- rbinom(k, sizeP, p)
m1 <- rbinom(k, sizeR, r)
dat <- cbind(n1=n1, n2=sizeP-n1, m1=m1, m2=sizeR-m1)
return(list( risk=p
, study=apply(dat, c(1,2), function(tmp) {max(1,tmp)})
))
}
#' Data generation: Sampling data of clinical trials
#'
#' Calculates log risk ratios from a study in the right format.
#'
#' @param study Study data of a clinical trial with
#' binomial outcomes.
#' @export
#' @examples
#' h_test <- .03
#' x_test <- cbind(1,1:13)
#' b_test <- c(0.02, 0.03)
#' s_test <- rep(2000, 13)
#' a_test <- rep(.3, 13)
#' rBinomGauss( h=h_test, s=s_test, a=a_test, r=0.03
#' , x=x_test, b=b_test)$study -> test
#' yvec(test)
#' dvec(test)
yvec <- function(study) {
return(log(study[,1]/(study[,1]+study[,2]))
- log(study[,3]/(study[,3]+study[,4])))
}
#' Data generation: Sampling data of clinical trials
#'
#' Calculates the variance estimate of log risk ratios from a study in
#' the right format. See the example below for details.
#'
#' @param study Study data of a clinical trial with
#' binomial outcomes.
#' @export
#' @examples
#' h_test <- .03
#' x_test <- cbind(1,1:13)
#' b_test <- c(0.02, 0.03)
#' s_test <- rep(2000, 13)
#' a_test <- rep(.3, 13)
#' rBinomGauss( h=h_test, s=s_test, a=a_test, r=0.03
#' , x=x_test, b=b_test)$study -> test
#' yvec(test)
#' dvec(test)
dvec <- function(study) {
return((1/study[,1]) - (1/(study[,1] + study[,2])) + (1/study[,3]) - (1/(study[,3]+ study[,4])))
}
#' Data generation: Gaussian-Gaussian model
#'
#' Random draws of response vectors y following the distribution
#' of a random effects meta regression model. Each column is
#' an independent draw.
#'
#' @param n number of draws.
#' @param h heterogeneity.
#' @param d heteroscedasticity.
#' @param x design matrix.
#' @param b regression coefficients.
#' @return A (k,n)-matrix. Each column is an independent draw.
#' @examples
#' x_test = cbind(1,1:13)
#' h_test = .03
#' d_test = rchisq(13, df=0.02)
#' b_test = c(0.02, 0.03)
#' rY(n=10, h=h_test, d=d_test, x=x_test, b=b_test)
#' @export
rY <- function( n # number of draws
, h # heterogeneity
, d # heteroscedasticity
, x # design matrix
, b # regression coefficients
) {
return(matrix(rnorm( n*length(d)
, mean=as.vector(x%*%b)
, sd=sqrt(h + d)), ncol=n))
}
#' Data generation: Gaussian-Gaussian model
#'
#' Random draws of heteroscedasticity responses of studies, where each
#' study in a random effects meta regression model follows a Gaussian
#' response. Thus D = (d * X) / (s-1) where X is chi-squared
#' distributed.
#'
#' @param n number of draws.
#' @param d heteroscedasticity.
#' @param s study sizes.
#' @return A (k,n)-matrix. Each column is an independent draw.
#' @examples
#' d_test = rchisq(13, df=0.02)
#' s_test = rep(100, 13)
#' rD(n=10, d=d_test, s=s_test)
#' @export
rD <- function( n # number of draws
, d # heteroscedasticity
, s # study sizes
) {
return(matrix(d * rchisq(n*length(s), df=(s-1)) / (s-1), ncol=n))
}
#' Data generation: Log-risk-ration of a binomial-Gaussian model
#'
#' Random draws of log risk ratios from a hierarchical binomial
#' Gaussian model.
#'
#' It is always assumed that at least one response in a study has
#' happend, i.e., a response of 0 in a treatment or control group is
#' rounded up to 1. Note that this may lead to an overestimation of
#' small risks. If possible, make sure your sample sizes are large
#' enough to compensate for this effect.
#'
#' @param n number of draws.
#' @param h heterogeneity.
#' @param s study sizes.
#' @param a balance of group assignments.
#' @param r fixed risk in the treatment group.
#' @param x design matrix.
#' @param b regression coefficients.
#' @return A (2k,n) matrix. Each column is an independent draw.
#' @examples
#' h_test <- .03
#' x_test <- cbind(1,1:13)
#' b_test <- c(0.02, 0.03)
#' s_test <- rep(2000, 13)
#' a_test <- rep(.3, 13)
#' rB(n=10, h=h_test, s=s_test, a=a_test, r=.3, x=x_test, b=b_test)
#' @export
rB <- function( n # number of draws
, h # heterogeneity
, s # study sizes
, a # balance of group assignments
, r # fixed risk in the treatment group
, x # design matrix
, b # regression coefficients
) {
func <- function() {
study <- rBinomGauss(h=h,s=s,a=a,r=r,x=x,b=b)$study
return(c(yvec(study=study), dvec(study=study)))
}
return(replicate(n, func()))
}
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.