# R/modelexamples.R In UdderQuarterInfectionData: Udder Quarter Infection Data

#### Documented in Gamma_Frailty_Interval_Censoring

globalVariables("udderquarterinfection")

#' @importFrom stats nlm
#' @importFrom utils data
NULL

#'@export
#'@title Gamma Frailty Interval Censoring
#'@description Application of the Gamma Frailty Interval Censoring Model on the Udder Quarter Infection Data Set. For more information see Details.
#'@param print.level Parameter of \code{\link[stats]{nlm}} (default=2): this argument determines the level of printing which is done during the minimization process. The default value of 0 means that no printing occurs, a value of 1 means that initial and final details are printed and a value of 2 means that full tracing information is printed.
#'
#'@details This function fits a parametric Weibull baseline hazard frailty model with gamma distributed frailties for the udder quarter infection data taking into consideration the interval censored nature of the data. Further theoretical details can be found in the paper in the reference
#'
#'@author Klara Goethals
#'@author Luc Duchateau
#'@references Goethals, K., Ampe, B., Berkvens, D., Laevens, H., Janssen, P. and Duchateau, L. (2009). Modeling interval-censored, clustered cow udder quarter infection times through the shared gamma frailty model. Journal of Agricultural Biological and Environmental Statistics 14, 1-14.
#'
#'@return Returns a list with the NLM result in \code{nlm} and the covariance matrix in \code{covmat}.
#'
#'@section R Code for Model :
#'The source R code for this model can found:
#'\itemize{
#'\item in the \code{doc/Models_R_Code.R} file in the package installation folder.
#'\item by accessing the function by calling \code{Gamma_Frailty_Interval_Censoring} (without brackets) or \code{getAnywhere("Gamma_Frailty_Interval_Censoring")}.
#'}
#'
#' @examples
#' \dontrun{
#' library(UdderQuarterInfectionData)
#' data("udderquarterinfection")
#'
#' Gamma_Frailty_Interval_Censoring()
#' # $nlm #' #$nlm$minimum #' # [1] 5670.491 #' # #' #$nlm$estimate #' # [1] 3.7967246 0.1201593 1.9672298 0.8590531 #' # #' #$nlm$gradient #' # [1] 0.0002924871 0.0017653292 -0.0005460029 0.0003265086 #' # #' #$nlm$hessian #' # [,1] [,2] [,3] [,4] #' # [1,] 23.22965 -117.7682 -39.93813 -10.10561 #' # [2,] -117.76825 15471.4753 567.24283 1228.87332 #' # [3,] -39.93813 567.2428 664.76359 24.63047 #' # [4,] -10.10561 1228.8733 24.63047 147.76479 #' # #' #$nlm$code #' # [1] 1 #' # #' #$nlm$iterations #' # [1] 22 #' # #' # #' #$covmat
#' # [,1]          [,2]          [,3]         [,4]
#' # [1,] 0.049281911  0.0001242730  0.0027853686  0.001872592
#' # [2,] 0.000124273  0.0001982213 -0.0001015391 -0.001623066
#' # [3,] 0.002785369 -0.0001015391  0.0017306214  0.000746460
#' # [4,] 0.001872592 -0.0016230660  0.0007464600  0.020269244
#' }
Gamma_Frailty_Interval_Censoring <- function(print.level=2){
data("udderquarterinfection",envir=environment())

upper=udderquarterinfection$right/91.31 # Divide by 91.31 to avoid small lambdas lower=udderquarterinfection$left/91.31
fail=udderquarterinfection$status cluster=udderquarterinfection$cowid
X=udderquarterinfection$lactation>1 #Calculate the number of clusters. clusternames <- levels(as.factor(cluster)) ncluster <- length(clusternames) # Create a udderquarterinfection with the variables cluster, the lower bound, the upper bound, # the censoring indicator and the covariate. udderquarterinfectionint <- as.matrix(cbind(cluster,lower,upper,fail,X)) # create subsets for right-censored and interval-censored observations cendata <- udderquarterinfectionint[udderquarterinfectionint[,4]==0,] intdata <- udderquarterinfectionint[udderquarterinfectionint[,4]==1,] # Create a list of signs that corresponds to the n_ik (here restricted to 4 events) signs <- list(1,c(1,-1)) for(i in 3:5) signs[[i]] <- kronecker(signs[[i-1]],c(1,-1)) # Function to calculate the loglikelihood per cluster CalcLogLikClust <- function(i,x) { theta <- x[1]; lambda <- x[2]; gamma <- x[3]; beta <- x[4] cenX <- cendata[cendata[,1]==clusternames[i],5] intL <- intdata[intdata[,1]==clusternames[i],2] nevents <- length(intL) crossprod <- 1 if(nevents>0){ #comment: if there are no events, crosspod=1 intX <- intdata[intdata[,1]==clusternames[i],5] # Calculate R*_ij and L*_ij intRster <- lambda*(intdata[intdata[,1]==clusternames[i],3]^gamma)*exp(intX*beta) intLster <- lambda*(intL^gamma)*exp(intX*beta) # Calculate the vector p_i crossprod <- c(exp(intLster[1]),exp(intRster[1])) if(nevents>1){ for(ik in 2:nevents) crossprod<-kronecker(crossprod,c(exp(intLster[ik]),exp(intRster[ik]))) } } # Loglikelihood for 1 cluster return( log(1/(theta^(1/theta))*sum((1/ ((sum(lambda*(as.vector(cendata[cendata[,1]==clusternames[i],3])^gamma) *exp(cenX*beta))+1/theta+log(crossprod))^(1/theta)))*signs[[nevents+1]])) ) } # Calculate full marginal loglikelihood (formula 5) CalcLogLik <- function(x) { -sum(sapply(1:ncluster,CalcLogLikClust,x=x)) } # Maximising the full marginal loglikelihood to obtain parameter estimates init <- c(1,1,1,1) results <- nlm(CalcLogLik,init,print.level=print.level, hessian=TRUE) # Can take a while! #$minimum
# [1] 5670.491
#
# $estimate # [1] 3.7967246 0.1201593 1.9672298 0.8590531 # #$gradient
# [1]  0.0002924871  0.0017653292 -0.0005460029  0.0003265086
#
# $hessian # [,1] [,2] [,3] [,4] # [1,] 23.22965 -117.7682 -39.93813 -10.10561 # [2,] -117.76825 15471.4753 567.24283 1228.87332 # [3,] -39.93813 567.2428 664.76359 24.63047 # [4,] -10.10561 1228.8733 24.63047 147.76479 # #$code
# [1] 1
#
# $iterations # [1] 22 # Calculate covariance matrix covmatr <- solve(results$hessian)
#             [,1]          [,2]          [,3]         [,4]
# [1,] 0.049281911  0.0001242730  0.0027853686  0.001872592
# [2,] 0.000124273  0.0001982213 -0.0001015391 -0.001623066
# [3,] 0.002785369 -0.0001015391  0.0017306214  0.000746460
# [4,] 0.001872592 -0.0016230660  0.0007464600  0.020269244

return(list(
nlm=results,
covmat=covmatr
))
}


## Try the UdderQuarterInfectionData package in your browser

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

UdderQuarterInfectionData documentation built on Sept. 6, 2017, 5:03 p.m.