# Extra functions: Water quality index ---------------------------------------- # Author: Ricardo R. Petterle/UFPR and Wagner H. Bonat LEG/UFPR --------------- # Create: 16/09/2017; Update: 14/10/2018 -------------------------------------- mc_non_aux <- function(n.resp) { position <- combn(n.resp, 2) list.Derivative <- list() n.par <- n.resp * (n.resp - 1)/2 for (i in 1:n.par) { Derivative <- matrix(0, ncol = n.resp, nrow = n.resp) Derivative[position[1, i], position[2, i]] <- Derivative[position[2, i], position[1, i]] <- 1 list.Derivative[i][[1]] <- Derivative } return(list.Derivative) } # Graph ----------------------------------------------------------------------- ac <- list(pad1=0.5, pad2=0.5, tck=0.5) mycol <- gray.colors(n=5) ps <- list(box.rectangle=list(col=1, fill=c("gray70")), box.umbrella=list(col=1, lty=1), dot.symbol=list(col=1), dot.line=list(col=1, lty=3), plot.symbol=list(col=1, cex=0.7), plot.line=list(col=1), plot.polygon=list(col="gray80"), superpose.line=list(col=mycol), superpose.symbol=list(col=mycol), par.main.text = list(font = 2, # make it bold just = "left", y = grid::unit(-2.9, "mm")), par.xlab.text = list(font = 1, # make it bold #just = "left", y = grid::unit(-0.2, "mm")), strip.background=list(col=c("gray90","gray70")), superpose.polygon=list(col=mycol), strip.background=list(col=c("gray90","gray70")), layout.widths=list( left.padding=0.25, right.padding=-1, ylab.axis.padding=0), layout.heights=list( bottom.padding=0.25, top.padding=0, axis.xlab.padding=0, xlab.top=0), axis.components=list(bottom=ac, top=ac, left=ac, right=ac) ) panel.segplotBy <- function(x, y, z, centers, subscripts, groups, f, ...){ d <- 2*((as.numeric(groups)-1)/(nlevels(groups)-1))-1 z <- as.numeric(z)+f*d panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...) } #' @title Delta Method #' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br} #' #' @description Compute standard errors for functions of model parameters #' using the delta method. #' #' @param g A function (string like formula) of model parameters. #' @param mean Vector of parameter estimates. #' @param cov Variance-covariance matrix. #' @param ses Logical. If TRUE returns the standard error, otherwise #' return the new variance-covariance. #' @keywords internal #' @details It is an internal function useful in general for summary #' function associated with Twin models. deltamethod <- function (g, mean, cov, ses = TRUE) { cov <- as.matrix(cov) n <- length(mean) if (!is.list(g)) g <- list(g) if ((dim(cov)[1] != n) || (dim(cov)[2] != n)) stop(paste("Covariances should be a ", n, " by ", n, " matrix")) syms <- paste("x", 1:n, sep = "") for (i in 1:n) assign(syms[i], mean[i]) gdashmu <- t(sapply(g, function(form) { as.numeric(attr(eval(deriv(form, syms)), "gradient")) })) new.covar <- gdashmu %*% cov %*% t(gdashmu) if (ses) { new.se <- sqrt(diag(new.covar)) new.se } else new.covar } # Odds ratio and CI ----------------------------------------------------------- mc_odds_ratio <- function(fit, names, response){ if(missing(response)) { response <- 1 } est <- coef(fit, type = "beta", std.error = T, response)$Estimates std <- coef(fit, type = "beta", std.error = T, response)$Std.error Z.values <- est/std p_value <- round((1-pnorm(abs(Z.values)))*2, dig = 3) Ic.min <- exp(est - qnorm(0.975) * std) Ic.max <- exp(est + qnorm(0.975) * std) out_temp <- data.frame("OR" = exp(est), "Ic.min" = Ic.min, "Ic.max" = Ic.max, "p-valor" = p_value) out <- round(out_temp[-1,], dig = 3) colnames(out) <- c("OR", "2.5%", "97.5%", "p-value") row.names(out) <- names return(out) } # Hat values ------------------------------------------------------------------ mc_hatvalues <- function(object, response, graph = FALSE){ if(missing(response)) { response <- 1 } # Setup X <- model.matrix(object$linear_pred[[response]][-2], object$data) N <- nrow(X) p <- ncol(X) N_resp <- length(object$linear_pred) N_respData <- dim(object$C)[1] ini = 1 fim = N_respData/N_resp if(response == 1) { ini <- 1 fim <- (N_respData/N_resp) } else if(response != 1) { ini <- response*(N_respData/N_resp) + 1 - (N_respData/N_resp) fim <- ini + (N_respData/N_resp) - 1 } # Covariance matrix (Weight matrix) W <- object$C[ini:fim, ini:fim] rootW <- matrix(0,N,N) l <- 1 while (l cut_point) if (length(extreme) > 0) for(i in 1:length(extreme)){ points(extreme[i], h[extreme[i]], pch = 20) text(extreme[i]+5, h[extreme[i]], extreme[i]) } } else return(h) } # Cook's distance ------------------------------------------------------------- mc_cookDist <- function(object, response, xlab, ylab, main, graph = FALSE){ if(missing(response)) { response <- 1 } # Setup X <- model.matrix(object$linear_pred[[response]][-2], object$data) N <- nrow(X) p <- ncol(X) N_resp <- length(object$linear_pred) N_respData <- dim(object$C)[1] mu <- fitted(object)[,response] y_observed <- object$observed[,response] error <- y_observed - mu ini = 1 fim = N_respData/N_resp if(response == 1) { ini <- 1 fim <- (N_respData/N_resp) } else if(response != 1) { ini <- response*(N_respData/N_resp) + 1 - (N_respData/N_resp) fim <- ini + (N_respData/N_resp) - 1 } # Covariance matrix (Weight matrix) W <- object$C[ini:fim, ini:fim] rootW <- matrix(0,N,N) l <- 1 while (l cut_point) if(length(extreme) > 0) for(i in 1:length(extreme)){ points(extreme[i], cdist[extreme[i]], pch = 20) text(extreme[i]+16, cdist[extreme[i]], extreme[i]) } } else return(cdist) } # DFBETAS --------------------------------------------------------------------- mc_DFBETA <- function(object, response, graph = FALSE){ if(missing(response)) { response <- 1 } # Setup X <- model.matrix(object$linear_pred[[response]][-2], object$data) N <- nrow(X) p <- ncol(X) N_resp <- length(object$linear_pred) N_respData <- dim(object$C)[1] mu <- fitted(object)[,response] y_observed <- object$observed[,response] error <- y_observed - mu ini = 1 fim = N_respData/N_resp if(response == 1) { ini <- 1 fim <- (N_respData/N_resp) } else if(response != 1) { ini <- response*(N_respData/N_resp) + 1 - (N_respData/N_resp) fim <- ini + (N_respData/N_resp) - 1 } # Hat values h <- mc_hatvalues(object, response, graph = FALSE) # Std error std.error <- sqrt(diag(object$C[ini:fim, ini:fim])) # DFBETAS XlX <- crossprod(X) iXlX <- solve(XlX) num <- t(sweep(iXlX%*%t(X), 2, error/(1 - h), "*")) den <- outer(std.error, sqrt(diag(iXlX)), "*") dfbetas <- num/den if(graph == TRUE){ par(mfrow = c(1, p)) for(i in 1:p){ plot(dfbetas[,i], main = colnames(dfbetas)[i], type = "h") } } else return(dfbetas) } # DFFITS ---------------------------------------------------------------------- mc_DFFITS <- function(object, response, graph = FALSE){ if(missing(response)) { response <- 1 } # Setup X <- model.matrix(object$linear_pred[[response]][-2], object$data) N <- nrow(X) p <- ncol(X) N_resp <- length(object$linear_pred) N_respData <- dim(object$C)[1] mu <- fitted(object)[,response] y_observed <- object$observed[,response] error <- y_observed - mu ini = 1 fim = N_respData/N_resp if(response == 1) { ini <- 1 fim <- (N_respData/N_resp) } else if(response != 1) { ini <- response*(N_respData/N_resp) + 1 - (N_respData/N_resp) fim <- ini + (N_respData/N_resp) - 1 } # Hat values h <- mc_hatvalues(object, response, graph = FALSE) # Std error std.error <- sqrt(diag(object$C[ini:fim, ini:fim])) # DFFITS s <- error/(std.error * sqrt((1 - h))) dffits <- s * sqrt((h/(1 - h))) if(graph == TRUE){ plot(dffits) } else return(dffits) } # Half-normal plot with simulated envelope ------------------------------------ mc_envelope <- function(object, response, envelope = "Hnormal", xlab, ylab, main, repl, pch){ if(missing(response)) { response <- 1 } # Setup linear_pred <- object$linear_pred[[response]] X <- model.matrix(linear_pred[-2], object$data) N <- nrow(X) p <- ncol(X) dataSet <- object$data N_resp <- length(object$linear_pred) N_respData <- dim(object$C)[1] ini = 1 fim = N_respData/N_resp if(response == 1) { ini <- 1 fim <- (N_respData/N_resp) } else if(response != 1) { ini <- response*(N_respData/N_resp) + 1 - (N_respData/N_resp) fim <- ini + (N_respData/N_resp) - 1 } # Stardardized residuals from fitted model ---------- # Covariance matrix (Weight matrix) W <- object$C[ini:fim, ini:fim] rootW <- matrix(0,N,N) l <- 1 while (l