# R/auxiliar_functions.R In iccCounts: Intraclass Correlation Coefficient for Count Data

#### Documented in ICCplot.GOFVarComp

#' @importFrom stats poisson pnbinom ppois predict qnorm runif simulate quantile
#' @importFrom VGAM pzipois pzinegbin
#' @importFrom gridExtra grid.arrange
#' @importFrom Deriv Deriv
#' @importFrom utils combn globalVariables
#' @import glmmTMB dplyr ggplot2

utils::globalVariables(c("m", "Diff", "lev1","lev2","Levels","per"))

# Transformations

f1<-function(s) exp(s)
f2<-function(s) exp(s)/(1+exp(s))
f3<-function(s) exp(s)^2
f4<-function(b1,b2) log(mu_f(b1,b2))

# Expectations and other parameters
mu_f<-function(b1,b2) f1(b1+f3(b2)/2)
mu_f2<-function(b1,b2,sb) f1(b1+(f3(b2)+sb)/2)

sa_f<-function(b2) f3(b2)
k_f<-function(b3) f1(b3) # gof parameter negbinom
p_f<-function(b4) f2(b4) # inflation parameter ZI

# Builds L matrix for 2x2 comparison of means
Lmat<-function(k){
aux1<-diag(k)*0
aux1[1,1:k]<-1

for (i in 1:(k-1)) aux1[1+i,1+i]=1

L=array(0,dim=c(k,k*(k-1)/2))
cont=0
for (i in 1:k){
for (j in 1:(k-1)){
if (i>j){
cont=cont+1
L[,cont]=aux1[,i]-aux1[,j]
}
}
}
return(L)
}

get_data<-function(data,y,id,met=NULL){

data<-as.data.frame(data)
vid<-,id]
vy<-,y]

if (is.null(met)==F){
vmet<-,met]
vmet<-as.factor(vmet)
x<-data.frame(vy,vid,vmet)
names(x)<-c("y","id","met")
if (length(unique(x\$met))<2) stop("Less than two methods")
return(x)
}  else {
x<-data.frame(vy,vid)
names(x)<-c("y","id")
return(x)
}

}

get_data_plot<-function(data,y,id,rm=NULL){

data<-as.data.frame(data)
vid<-,id]
vy<-,y]

if (is.null(rm)==F){
vrm<-,rm]
vrm<-as.factor(vrm)
x<-data.frame(vy,vid,vrm)
names(x)<-c("y","id","rm")
if (length(unique(x\$rm))<2) stop("Less than two repeated measures")
return(x)
}  else {
x<-data.frame(vy,vid)
names(x)<-c("y","id")
return(x)
}

}

sb_f<-function(k,b,varB){

L<-Lmat(k)   # Building L matrix

difmed=t(L)%*%b  # Calculating observers sum of squares

A<-L%*%t(L)
aux2<-(t(difmed)%*%difmed)-sum(diag(A%*%varB))
sb<-max(aux2/(k*(k-1)),0)
return(sb)
}

# Variance of sb
var_sb<-function(k,b,varB){
L<-Lmat(k)    # Building L matrix
A<-L%*%t(L)
xx<-( 2*sum(diag( ( (A%*%varB)**2 ) ) )+ 4*( t(b)%*%A%*%varB%*%A%*%b ) ) /((k*(k-1))**2)
return(xx)
}

ic_r<-function(r,var.r,conf=0.95){

alpha<-1-conf
z<-0.5*log((1+r)/(1-r))
var.z<-var.r/(((1+r)**2)*((1-r)**2))
ic.z<-z+c(-1,1)*as.vector(qnorm(1-(alpha/2))*sqrt(var.z))
ic.r<-(exp(2*ic.z)-1)/(exp(2*ic.z)+1)
out<-matrix(c(r,sqrt(var.r),ic.r),nrow=1)

colnames(out)<-c("ICC","SE ICC",paste(conf*100,"% CI LL",sep=""),
paste(conf*100,"% CI UL",sep=""))
return(out)

}

varB<-function(model,k) model\$sdr\$cov.fixed[1:k,1:k]

#' Prints the intraclass correlation coefficient
#'
#' Prints the intraclass correlation coefficient (ICC)
#'
#' @param x An object of class *iccc*
#' @param digits Number of digits to print
#' @return A vector with the ICC estimate, its standard error and confidence interval.
#' @export
#' @examples
#' # Poisson model. Repeatability setting.
#' iccpois<-icc_counts(Grimso,y="Tot",id="TransectID")
#' ICC(iccpois)
#' # Negative Binomial with proportional extra-dispersion. Concordance setting
#' iccnb2<-icc_counts(AF,y="y",id="id",met="met",type="con",fam="nbinom2")
#' ICC(iccnb2)
#' # Zero-inflated Poisson model. Repeatability setting
#' icczip<-icc_counts(EPP,y="Social",id="id",fam="zip")
#' ICC(icczip)
#' @seealso [icc_counts()]
ICC<-function(x,digits = getOption("digits")){

fam<-x\$model\$modelInfo\$family\$family
zi<-x\$model\$modelInfo\$allForm\$ziformula
if (zi==~1) ot<-"zero inflated" else ot<-""

message("Model: ",paste(fam,ot))

return(round(x\$ICC,digits = digits))

}

#' GLMM variance components
#'
#' Prints the GLMM variance components and related parameters to estimate the intraclass correlation coefficient (ICC)
#'
#' @param x An object of class *iccc*
#' @param digits Number of digits to print
#' @return A vector with the variance components and related parameters
#' @export
#' @examples
#' # Poisson model. Repeatability setting.
#' iccpois<-icc_counts(Grimso,y="Tot",id="TransectID")
#' VarComp(iccpois)
#' # Negative Binomial with proportional extra-dispersion. Concordance setting
#' iccnb2<-icc_counts(AF,y="y",id="id",met="met",type="con",fam="nbinom2")
#' VarComp(iccnb2)
#' # Zero-inflated Poisson model. Repeatability setting
#' icczip<-icc_counts(EPP,y="Social",id="id",fam="zip")
#' VarComp(icczip)
#' @seealso [icc_counts()]
VarComp<-function(x,digits = getOption("digits")){

fam<-x\$model\$modelInfo\$family\$family
zi<-x\$model\$modelInfo\$allForm\$ziformula
if (zi==~1) ot<-"zero inflated" else ot<-""
message("Model: ",paste(fam,ot))
return(round(x\$varcomp,digits = digits))

}

count_zero<-function(x) sum(x==0)

#' Goodness of fit plots
#'
#' Draws the plots to assess the goodness of fit
#'
#' @param x An object of class *GOF* generated by *GOF_check* function.
#' @param type Which plot to draw. Values: *all* (default); *envelope* for envelopes of randomized quantile residuals; *dispersion* for plot to assess the dispersion; *zeros* for plot to assess the zero inflation.
#' @param ... Ignore
#' @method plot GOF
#' @export
#' @examples
#' \donttest{
#' # Poisson model. Repeatability setting.
#' iccpois<-icc_counts(EPP,y="Social",id="id")
#' iccpois.gof<-GOF_check(iccpois)
#' plot(iccpois.gof)
#' plot(iccpois.gof,type="envelope")
#' plot(iccpois.gof,type="dispersion")
#' plot(iccpois.gof,type="zeros")
#' }
#' @seealso [GOF_check()]
plot.GOF<-function(x,type=c("all","envelope","dispersion","zeros"),...){
type<-match.arg(type,choices=c("all","envelope","dispersion","zeros"))

if (type=="all") grid.arrange(x\$plot_env,x\$plot_var,x\$plot_zi)
if (type=="envelope") return(x\$plot_env)
if (type=="dispersion") return(x\$plot_var)
if (type=="zeros") return(x\$plot_zi)

}

## Try the iccCounts package in your browser

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

iccCounts documentation built on May 29, 2024, 7:57 a.m.