# Unit gamma mixed models for continuous bounded data ------------------------- # Analysis: Water quality index data ------------------------------------------ # Author: Ricardo Rasmussen Petterle/UFPR ------------------------------------- # Date: 05/09/2020 ------------------------------------------------------------ rm(list=ls()) # BEGIN HEADER ---------------------------------------------------------------- # Loading extra packages require(plyr) require(TMB) # Loading data set dados <- read.table("data_set2_UG.txt", h = T) # Preparing data set dados$TRIM <- as.factor(dados$TRIM) levels(dados$TRIM) dados$LOCAL <- factor(dados$LOCAL, levels = levels(dados$LOCAL)[c(2,3,1)]) levels(dados$LOCAL) <- c("Upstream", "Reservoir", "Downstream") dados$id <- factor(dados$id) levels(dados$id) <- paste("U", 1:16, sep = "") dados <- na.exclude(dados) # END HEADER ------------------------------------------------------------------ # Fitting # Model 1: unit gamma with random intercept ----------------------------------- compile("unit_gamma_fixed.cpp") dyn.load(dynlib("unit_gamma_fixed")) # Data X <- model.matrix(~ TRIM + LOCAL, dados) X <- as(X,"dgTMatrix") Z <- model.matrix(~ factor(id)-1, dados) Z <- as(Z,"dgTMatrix") betas <- rnorm(ncol(X)) u <- rnorm(ncol(Z)) data <- list(y = dados$y, X = X, Z = Z) ## Initial values parameters <- list(beta = betas*0, logphi = 1, logsdu = 1, u = u*0) # Log-likelihood function obj <- MakeADFun(data, parameters, DLL = "unit_gamma_fixed", random = "u", hessian = TRUE, silent = TRUE) # Optimization model1 <- nlminb(obj$par, obj$fn, obj$gr) temp <- sdreport(obj) summary(temp, "fixed", p.value = T) summary(temp, "report", p.value = T) # Model 2: beta with random intercept ----------------------------------------- compile("beta_fixed.cpp") dyn.load(dynlib("beta_fixed")) # Log-likelihood function obj <- MakeADFun(data, parameters, DLL = "beta_fixed", random = "u", hessian = TRUE, silent = TRUE) # Optimization model2 <- nlminb(obj$par, obj$fn, obj$gr) temp <- sdreport(obj) summary(temp, "fixed", p.value = T) summary(temp, "report", p.value = T) # Model 3: unit gamma nested model -------------------------------------------- compile("unit_gamma_nested.cpp") dyn.load(dynlib("unit_gamma_nested")) # Data X <- model.matrix(~ TRIM + LOCAL, dados) X <- as(X,"dgTMatrix") Z <- model.matrix(~ factor(id)-1, dados) Z <- as(Z,"dgTMatrix") W <- model.matrix(~ TRIM:factor(id)-1, dados) W <- as(W,"dgTMatrix") betas <- rnorm(ncol(X)) ui <- rnorm(ncol(Z)) uij <- rnorm(ncol(W)) data <- list(y = dados$y, X = X, Z = Z, W = W) ## Initial values parameters <- list(beta = betas*0, logphi = 1, logtau1 = 1, logtau2 = 1, ui = ui*0, uij = uij*0) # Log-likelihood function obj <- MakeADFun(data, parameters, DLL = "unit_gamma_nested", random = c("ui","uij"), hessian = TRUE, silent = TRUE) # Optimization model3 <- nlminb(obj$par, obj$fn, obj$gr) temp <- sdreport(obj) summary(temp, "fixed", p.value = TRUE) summary(temp, "report", p.value = TRUE) # Model 4: beta nested model -------------------------------------------------- compile("beta_nested.cpp") dyn.load(dynlib("beta_nested")) # Log-likelihood function obj <- MakeADFun(data, parameters, DLL = "beta_nested", random = c("ui","uij"), hessian = TRUE, silent = TRUE) # Optimization model4 <- nlminb(obj$par, obj$fn, obj$gr) temp <- sdreport(obj) summary(temp, "fixed", p.value = TRUE) summary(temp, "report", p.value = TRUE) # END -------------------------------------------------------------------------