R/single.r

Defines functions singletable

Documented in singletable

######################################################################################
### Purpose: Create object "singletable"
### Input:   data(y1,n1,y2,n2), hyperparemeters(a1,a2,b1,b2,rho),
###          comaparative measure, model, method, significance level (alpha),
###          number of samples (nsam)
### Output:  S3 object "singletable", refer to the help file of singletable
### Author:  Sheng Luo, Yong Chen, Xiao Su, Haitao Chu
### Data:    7/13/2012
######################################################################################
singletable <- function(y1=y1,n1=n1,y2=y2,n2=n2,measure=measure,model="Sarmanov",
                        method="exact",a1=0.5,b1=0.5,a2=0.5,b2=0.5,rho=0,alpha=0.05,
                        nsam=10000) {
  
  if (measure=="RD"& method=="exact") {
    print("only sampling based mehtod is available for RD")
    method <- "sampling"
  }
  if (!model%in%c("Sarmanov","Independent")) stop("Only \'Independent\' and \'Sarmanov\' are supported. ")

  if(length(y1)>=2 |length(n1)>=2 |length(y2)>=2 |length(n2)>=2 )
    stop ("only for single table analysis \n")

  if(y1<0|y2<0) stop("y1,y2 should be greater than 0")
  if(n1<=0|n2<=0) stop("n1,n2 should not be less or equal to 0")
   
  if ((trunc(y1)!=y1)|(trunc(n1)!=n1)|(trunc(y2)!=y2)|(trunc(n2)!=n2))
    stop("y1,n1,y2,n2 should be interger")
  
  if(y1>n1) stop("y1 should be less than n1")
  if(y2>n2) stop("y2 should be less than n2")

  if (any(y1==n1)) {
    index=which(y1==n1)
    for(i in 1:length(index)) {
      y1[index[i]]=y1[index[i]]-0.02
    }
  }

  if (any(y2==n2)) {##if y1=n2, the bugs will down
    index=which(y2==n2)
    for(i in 1:length(index)) {
      y2[index[i]]=y2[index[i]]-0.02
    }
  }	   
   
  # check the range of parameters
  if(model=="Sarmanov") {
    cc <- sqrt(a1*a2*b1*b2)/sqrt((a1+b1+1)*(a2+b2+1))
    upper.bound <- cc/max(a1*b2, a2*b1)
    lower.bound <- -cc/max(a1*a2, b1*b2)
    rho.range<-c(lower.bound,upper.bound)
    names(rho.range)<-c("lower.bound","upper.bound")
    cat("Range of corelation (rho)",fill=T)
    cat("Lower bound:", rho.range[1],fill=T)
    cat("Upper bound:",rho.range[2],fill=T)
    if (rho > upper.bound | rho < lower.bound) stop(paste("rho is out of bound: ",
    lower.bound, upper.bound))
  }
  
  #save parameter
  if(model=="Independent") {
    parameter <- c(a1,b1,a2,b2)
    names(parameter) <- c("a1","b1","a2","b2")
  }

  if(model=="Sarmanov") {
    parameter<-c(a1,b1,a2,b2,rho)
    names(parameter) <- c("a1","b1","a2","b2","rho")
  }
                 
  # save data
  dataset <- c(y1,n1,y2,n2)
  names(dataset) <- c("y1","n1","y2","n2")

  if (measure=="OR") measurename<-"Odds ratio"
  if (measure=="RR") measurename<-"Relative risk"
  if (measure=="RD") measurename<-"Risk difference"

  mysample <- list()
  temp <- sampling(a1=a1,b1=b1,a2=a2,b2=b2,rho=rho,n1=n1,y1=y1,
                            n2=n2,y2=y2,measure=measure,model=model,nsam=nsam)
  upper=quantile(temp,prob=0.9995,na.rm=T)
  mysample[[1]]=temp[temp<=upper]  
  temp<- sampling(a1=a1,b1=b1,a2=a2,b2=b2,rho=rho,n1=0,y1=0,n2=0,
                            y2=0,measure=measure,model=model,nsam=nsam)
  upper=quantile(temp,prob=0.9995,na.rm=T)
  mysample[[2]]=temp[temp<=upper]
  
  dens <- list()
  if (method=="sampling") {
  #emprial density
    if (measure!="RD") dens <- lapply(mysample,density,from=0,n=2048)
    if (measure=="RD") dens <- lapply(mysample,density,from=-1,to=1,n=2048)
  }

  if (method=="exact") {
    xmin.post <- quantile(mysample[[1]],prob=0.0025,na.rm=T)
    xmax.post <- quantile(mysample[[1]],prob=0.9975,na.rm=T)
    dens[[1]] <- dens.post(y1=y1, n1=n1, y2=y2, n2=n2, a1=a1, b1=b1, a2=a2, b2=b2,
                           rho=rho,grid.start=xmin.post,grid.end=xmax.post,grid.num=1000,
                           measure=measure,model=model)
    dens[[2]] <- dens.post(y1=0, n1=0, y2=0, n2=0, a1=a1, b1=b1, a2=a2, b2=b2,
                           rho=rho, grid.start=xmin.post, grid.end=xmax.post, grid.num=1000,
                           measure=measure,model=model)
  }
  names(dens) <- c("Posterior","Prior")
  studynames <- c("Posterior","Prior")
  result <- list(measure=measure,model=model,method=method,dataset=dataset,
                 parameter=parameter,alpha=alpha,density=dens,sample=mysample,
                 studynames=studynames,measurename=measurename)

  class(result) <- "singletable"
  invisible(result)
}     

Try the mmeta package in your browser

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

mmeta documentation built on May 29, 2017, 8:28 p.m.