# Data analyses: Flexible regression models for continuous bounded data # Data analysis 1 - IQA data ----------------------------------------- # Author: Wagner Hugo Bonat LEG/UFPR --------------------------------- # Date: 26/10/2017 --------------------------------------------------- rm(list=ls()) # Loading extra packages require(mcglm) require(Matrix) require(betareg) require(simplexreg) require(VGAM) require(pracma) # Loading extra functions source("functions.R") # Data set 1 - Water quality index ----------------------------------- data1 <- read.table("../DataSets/IQAData.txt", header = TRUE) data1$TRIM <- as.factor(data1$TRIM) data1$LOCAL <- factor(data1$LOCAL, levels=c("MONT","RESER","JUSA")) levels(data1$LOCAL) <- c("Upstream","Reservoir","Downstream") par(mfrow = c(1,3), mar=c(2.6, 2.5, 1, 0.1), mgp = c(1.6, 0.6, 0)) hist(data1$IQA, ylab = "Frequency", xlab = "IQA", main = "A") boxplot(data1$IQA ~ data1$LOCAL, ylab = "IQA", xlab = "LOCAL", notch = TRUE, main = "B") boxplot(data1$IQA ~ data1$TRIM, ylab = "IQA", xlab = "QUARTER", notch = TRUE, main = "C") # Linear predictor form1 <- IQA ~ LOCAL + TRIM Z0_1 <- mc_id(data1) # Beta fit ----------------------------------------------------------- fit_beta1 <- betareg(form1, data = data1) # Simplex fit -------------------------------------------------------- fit_sim1 <- vglm(form1, family = "simplex", data = data1) # McGLM fit ---------------------------------------------------------- fit_mc1 <- mcglm(c(form1), list(Z0_1), link = "logit", variance = "binomialP", power_fixed = FALSE, control_algorithm = list(tunning = 0.75, correct = FALSE), data = data1) logLik(fit_beta1) logLik(fit_sim1) # Pseudo-logLik beta fit pll_beta <- plogLik_iid(pred = fitted(fit_beta1), y = data1$IQA, dispersion = coef(fit_beta1)[7], model = "beta") # Pseudo-logLik simplex fit pll_sim <- plogLik_iid(pred = fitted(fit_sim1), y = data1$IQA, dispersion = exp(coef(fit_sim1)[2]), model = "simplex") # Pseudo-logLik mcglm fit plogLik(fit_mc1) # Estimates est_beta <- coef(fit_beta1)[-7] est_sim <- coef(fit_sim1)[-2] est_mc <- coef(fit_mc1)$Estimates[-c(7,8)] std_beta <- sqrt(diag(vcov(fit_beta1)[-7,-7])) std_sim <- sqrt(diag(vcov(fit_sim1)[-2,-2])) std_mc <- coef(fit_mc1, std.error = TRUE)$Std.error[-c(7,8)] z_beta <- est_beta/std_beta z_sim <- est_sim/std_sim z_mc <- est_mc/std_mc # Comparison Table round(cbind(est_beta, std_beta, z_beta, est_sim, std_sim, z_sim, est_mc, std_mc, z_mc),2) c(logLik(fit_beta1),logLik(fit_sim1)) c(pll_beta, pll_sim, plogLik(fit_mc1)$plogLik) summary(fit_mc1) summary(fit_beta1) mc_conditional_test(fit_mc1, parameters = c("power11","tau11"), test = 1, fixed = 2) mc_conditional_test(fit_mc1, parameters = c("power11","tau11"), test = 2, fixed = 1)