######################################################################## # Final bivariate count model - Hunter data ############################ # Author: Wagner Hugo Bonat LEG/IMADA ################################## # Date:12/08/2015 ###################################################### ######################################################################## rm(list=ls()) # Loading extra packages ----------------------------------------------- require(mcglm) require(Matrix) require(plyr) require(lattice) require(latticeExtra) require(splines) # Loading data set ----------------------------------------------------- setwd("~/Dropbox/Doutorado/HuntingBioko2016/DataAnalysis") data <- read.table("HuntingData120816.txt", header = TRUE) # Factor data$Alt <- as.factor(data$Alt) data$Sex <- as.factor(data$Sex) data$Method <- as.factor(data$Method) # Hunter day effect data$OBS <- 1:dim(data)[1] data$HunterDay <- as.factor(paste(data$Hunter, data$Month, sep = "")) data$HunterDay <- droplevels(data$HunterDay) data <- arrange(data, Hunter) # Number of observations per Hunter range(table(data$Hunter)) # Number of observation per Hunter and month range(table(data$HunterDay)) # Matrix linear predictor ---------------------------------------------- Z0 <- Diagonal(dim(data)[1], 1) mp1.ot <- list(Z0) # Selecting the linear predictor for OT -------------------------------- # Iteration 1 ---------------------------------------------------------- form1.ot <- OT ~ 1 fit1.ot <- mcglm(linear_pred = c(form1.ot), matrix_pred = list(mp1.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(TRUE), offset = list(log(data$Offset)), data = data) summary(fit1.ot) scope <- c("Sex","Method","Alt","poly(Month,4)", "Sex*Method","Sex*Alt", "Sex*poly(Month,4)","Method*Alt", "Method*poly(Month,4)", "Alt*poly(Month,4)") SIC1 <- mc_sic(fit1.ot, scope = scope, data=data, response = 1, penalty = 2) SIC1 # Iteration 2 ---------------------------------------------------------- form2.ot <- OT ~ Method*Alt fit2.ot <- mcglm(linear_pred = c(form2.ot), matrix_pred = list(mp1.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit2.ot) anova(fit2.ot) scope <- c("Sex","poly(Month,4)", "Sex*Method","Sex*Alt", "Sex*poly(Month,4)","Method*poly(Month,4)", "Alt*poly(Month,4)") SIC2 <- mc_sic(fit2.ot, scope = scope, data=data, response = 1, penalty = 2) SIC2 # Iteration 3 ---------------------------------------------------------- form3.ot <- OT ~ Method*Alt + Sex*Alt fit3.ot <- mcglm(linear_pred = c(form3.ot), matrix_pred = list(mp1.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit3.ot) anova(fit3.ot) scope <- c("poly(Month,4)", "Sex*Method", "Sex*poly(Month,4)","Method*poly(Month,4)", "Alt*poly(Month,4)") SIC3 <- mc_sic(fit3.ot, scope = scope, data=data, response = 1, penalty = 2) SIC3 # Iteration 4 ---------------------------------------------------------- form4.ot <- OT ~ Method*Alt + Sex*Alt + Alt*poly(Month,4) fit4.ot <- mcglm(linear_pred = c(form4.ot), matrix_pred = list(mp1.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit4.ot) anova(fit4.ot) scope <- c("Sex*Method", "Sex*poly(Month,4)","Method*poly(Month,4)") SIC4 <- mc_sic(fit4.ot, scope = scope, data=data, response = 1, penalty = 2) SIC4 # Final model OT -------------------------------------------------------- final.ot <- OT ~ Method*Alt + Sex*Alt + Alt*poly(Month, 4) fit.ot <- mcglm(linear_pred = c(final.ot), matrix_pred = list(mp1.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit.ot) anova(fit.ot) # END ------------------------------------------------------------------ # Selecting the linear predictor for BD -------------------------------- # Iteration 1 ---------------------------------------------------------- mp1.bd <- list(Z0) form1.bd <- BD ~ 1 fit1.bd <- mcglm(linear_pred = c(form1.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(TRUE), offset = list(log(data$Offset)), data = data) summary(fit1.bd) scope <- c("Sex","Method","Alt","poly(Month,3)", "Sex*Method","Sex*Alt","Sex*poly(Month,3)", "Method*Alt","Method*poly(Month,3)","Alt*poly(Month,3)") mc_sic(fit1.bd, scope = scope, data=data, response = 1, penalty = 2) # Iteration 2 ---------------------------------------------------------- form2.bd <- BD ~ Method*Alt fit2.bd <- mcglm(linear_pred = c(form2.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit2.bd) anova(fit2.bd) scope <- c("Sex","poly(Month,3)", "Sex*Method","Sex*Alt","Sex*poly(Month,3)", "Method*poly(Month,3)","Alt*poly(Month,3)") mc_sic(fit2.bd, scope = scope, data=data, response = 1, penalty = 2) # Iteration 3 ---------------------------------------------------------- form3.bd <- BD ~ Method*Alt + Sex*Alt fit3.bd <- mcglm(linear_pred = c(form3.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit3.bd) anova(fit3.bd) form31.bd <- BD ~ Method*Alt + Sex fit31.bd <- mcglm(linear_pred = c(form31.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit31.bd) anova(fit31.bd) scope <- c("poly(Month,3)", "Sex*Method","Sex*poly(Month,3)", "Method*poly(Month,3)","Alt*poly(Month,3)") mc_sic(fit31.bd, scope = scope, data=data, response = 1, penalty = 2) # Iteration 4 ---------------------------------------------------------- form4.bd <- BD ~ Method*Alt + Sex + Alt*poly(Month, 3) fit4.bd <- mcglm(linear_pred = c(form4.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit4.bd) anova(fit4.bd) scope <- c("Sex*Method","Sex*poly(Month,3)", "Method*poly(Month,3)") mc_sic(fit4.bd, scope = scope, data=data, response = 1, penalty = 2) # Iteration 5 ---------------------------------------------------------- form5.bd <- BD ~ Method*Alt + Sex + Alt*poly(Month, 3) + Sex*Method fit5.bd <- mcglm(linear_pred = c(form5.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit5.bd) anova(fit5.bd) scope <- c("Sex*poly(Month,3)", "Method*poly(Month,3)") mc_sic(fit5.bd, scope = scope, data=data, response = 1, penalty = 2) anova(fit5.bd) # Final linear predictor for BD ---------------------------------------- final.bd <- BD ~ Method*Alt + Sex + Alt*poly(Month, 3) fit.bd <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) summary(fit.bd) anova(fit.bd) #----------------------------------------------------------------------- # Components for the matrix linear predictor --------------------------- data.hunter <- split(data, data$Hunter) # Known matrices ------------------------------------------------------- Z0 <- Diagonal(dim(data)[1],1) Z_hunter <- list() Z_long <- list() Z_HM <- list() for(i in 1:52) { Z <- model.matrix(~ Alt + Method + Sex, data = data.hunter[[i]]) Z_temp <- list() for(j in 1:dim(Z)[2]) { Z_temp[[j]] <- Z[,j]%*%t(Z[,j]) } Z_hunter[[i]] <- Z_temp temp <- 1/as.matrix(dist(data.hunter[[i]]$Month, diag = TRUE, upper = TRUE)) temp[temp == Inf] <- 0 temp <- Matrix(temp, sparse = TRUE) Z_long[[i]] <- temp tt <- table(data.hunter[[i]]$HunterDay) Z_HM_temp <- list() for(k in 1:length(tt)){ Z_HM_temp[[k]] <- tcrossprod(rep(1, tt[k])) } Z_HM[[i]] <- bdiag(Z_HM_temp) } Z_list <- list() for(i in 1:7) { Z_list[[i]] <- bdiag(lapply(Z_hunter, function(x,i)x[[i]], i = i)) } image(Z_list[[1]]) ZHM <- bdiag(Z_HM) image(ZHM) image(Z_list[[5]]) # Longitudinal effect -------------------------------------------------- Zlong <- bdiag(Z_long) # Selecting the matrix linear predictor for OT ------------------------- # Iteration 1 ---------------------------------------------------------- mp1.ot <- list(Z0) fit.ot1 <- mcglm(linear_pred = c(final.ot), matrix_pred = list(mp1.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) anova(fit.ot1) scope <- c(Z_list, ZHM, Zlong) id <- 1:9 mc_sic_covariance(fit.ot1, scope = scope, idx = id, data=data, response=1) # Iteration 2 ---------------------------------------------------------- mp2.ot <- list(Z0, Zlong) fit.ot2 <- mcglm(linear_pred = c(final.ot), matrix_pred = list(mp2.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), control_algorithm = list(max_iter = 100, tunning = 0.5), data = data) anova(fit.ot2) summary(fit.ot2) scope <- c(Z_list, ZHM) id <- 1:8 mc_sic_covariance(fit.ot2, scope = scope, idx = id, data=data, response=1) # Iteration 3 ---------------------------------------------------------- mp3.ot <- list(Z0, Zlong, ZHM) fit.ot3 <- mcglm(linear_pred = c(final.ot), matrix_pred = list(mp3.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), control_algorithm = list(max_iter = 100, tunning = 0.5), data = data) anova(fit.ot3) summary(fit.ot3) scope <- c(Z_list) id <- 1:7 mc_sic_covariance(fit.ot3, scope = scope, idx = id, data=data, response=1) # Iteration 4 ---------------------------------------------------------- mp4.ot <- list(Z0, ZHM, Z_list[[1]]) fit.ot4 <- mcglm(linear_pred = c(final.ot), matrix_pred = list(mp4.ot), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), control_algorithm = list(max_iter = 100, tunning = 0.5), data = data) anova(fit.ot4) summary(fit.ot4) scope <- c(Z_list[[2]],Z_list[[3]],Z_list[[4]],Z_list[[5]],Z_list[[6]],Z_list[[7]]) id <- 1:6 mc_sic_covariance(fit.ot4, scope = scope, idx = id, data=data, response=1) # Final model OT ------------------------------------------------------- mp.ot.final <- list(Z0, ZHM) fit.ot.final <- mcglm(linear_pred = c(final.ot), matrix_pred = list(mp.ot.final), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), control_algorithm = list(max_iter = 100, tunning = 1), data = data) summary(fit.ot.final) anova(fit.ot.final) #----------------------------------------------------------------------- # Selecting the matrix linear predictor for BD ------------------------- # Iteration 1 ---------------------------------------------------------- mp1.bd <- list(Z0) fit.bd1 <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp1.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) anova(fit.bd1) scope <- c(Z_list, ZHM, Zlong) id <- 1:9 mc_sic_covariance(fit.bd1, scope = scope, idx = id, data=data, response=1) # Iteration 2 ---------------------------------------------------------- mp2.bd <- list(Z0, ZHM) fit.bd2 <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp2.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) anova(fit.bd2) summary(fit.bd2) scope <- c(Z_list, Zlong) id <- 1:8 mc_sic_covariance(fit.bd2, scope = scope, idx = id, data=data, response=1) # Iteration 3 ---------------------------------------------------------- mp3.bd <- list(Z0, ZHM, Z_list[[6]]) fit.bd3 <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp3.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) anova(fit.bd3) summary(fit.bd3) scope <- c(Z_list[[1]], Z_list[[2]],Z_list[[3]],Z_list[[4]],Z_list[[5]], Z_list[[7]], Zlong) id <- 1:7 mc_sic_covariance(fit.bd3, scope = scope, idx = id, data=data, response=1) # Iteration 4 ---------------------------------------------------------- mp4.bd <- list(Z0, ZHM, Z_list[[6]], Z_list[[2]]) fit.bd4 <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp4.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) anova(fit.bd4) summary(fit.bd4) scope <- c(Z_list[[1]], Z_list[[3]],Z_list[[4]],Z_list[[5]], Z_list[[7]], Zlong) id <- 1:6 mc_sic_covariance(fit.bd4, scope = scope, idx = id, data=data, response=1) # Iteration 5 ---------------------------------------------------------- mp5.bd <- list(Z0, ZHM, Z_list[[6]], Zlong) fit.bd5 <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp5.bd), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), data = data) anova(fit.bd5) summary(fit.bd5) scope <- c(Z_list[[1]], Z_list[[3]],Z_list[[4]],Z_list[[5]], Z_list[[7]]) id <- 1:5 mc_sic_covariance(fit.bd5, scope = scope, idx = id, data=data, response=1) # END mp.bd.final <- list(Z0, ZHM, Z_list[[6]], Zlong) fit.bd.final <- mcglm(linear_pred = c(final.bd), matrix_pred = list(mp.bd.final), link = c("log"), variance = c("poisson_tweedie"), power_fixed = c(FALSE), offset = list(log(data$Offset)), control_algorithm = list(max_iter = 30), data = data) anova(fit.bd.final) summary(fit.bd.final) #----------------------------------------------------------------------- # Fitting the joint bivariate model ------------------------------------ fit.jt <- mcglm(c(final.bd, final.ot), list(mp.bd.final, mp.ot.final), link = c("log","log"), variance = c("poisson_tweedie","poisson_tweedie"), power_fixed = c(FALSE, FALSE), offset = list(log(data$Offset), log(data$Offset)), data = data, control_algorithm = list(tunning = 0.5, verbose = FALSE, max_iter = 100)) summary(fit.jt, print = c("power","Dispersion")) anova(fit.jt) # Additional models - bs splines instead of poly ----------------------- require(splines) final.bd.bs <- BD ~ Method*Alt + Sex + Alt*bs(Month, 3) final.ot.bs <- OT ~ Method*Alt + Sex*Alt + Alt*bs(Month, 4) fit.jt.bs <- mcglm(c(final.bd.bs, final.ot.bs), list(mp.bd.final, mp.ot.final), link = c("log","log"), variance = c("poisson_tweedie","poisson_tweedie"), power_fixed = c(FALSE, FALSE), offset = list(log(data$Offset), log(data$Offset)), data = data, control_algorithm = list(tunning = 0.5, tol = 1e-04, verbose = FALSE, max_iter = 500)) anova(fit.jt.bs) summary(fit.jt.bs) # Final model using only linear effect of Month effect final.bd.linear <- BD ~ Method*Alt + Sex + Alt*Month final.ot.linear <- OT ~ Method*Alt + Sex*Alt + Alt*Month fit.jt.linear <- mcglm(c(final.bd.linear, final.ot.linear), list(mp.bd.final, mp.ot.final), link = c("log","log"), variance = c("poisson_tweedie","poisson_tweedie"), power_fixed = c(FALSE, FALSE), offset = list(log(data$Offset), log(data$Offset)), data = data, control_algorithm = list(tunning = 0.5, tol = 1e-02, verbose = FALSE, max_iter = 500)) anova(fit.jt.linear) summary(fit.jt.linear) # Quadratic poly for BD ------------------------------------------------ mp.bd.final1 <- list(Z0, ZHM, Zlong) final.bd.quad <- BD ~ Method*Alt + Sex + Alt*poly(Month,2) final.ot.quad <- OT ~ Method*Alt + Sex*Alt + Alt*poly(Month,4) fit.jt.quad <- mcglm(c(final.bd.quad, final.ot.quad), list(mp.bd.final1, mp.ot.final), link = c("log","log"), variance = c("poisson_tweedie","poisson_tweedie"), power_fixed = c(FALSE, FALSE), offset = list(log(data$Offset), log(data$Offset)), data = data, control_algorithm = list(tunning = 0.25, method = "chaser", tol = 1e-04, verbose = TRUE, max_iter = 500)) # Comparing the additional models # Model presented in the paper gof(fit.jt) # Model using bs splines gof(fit.jt.bs) # Model using only linear trend gof(fit.jt.linear) # Model using quadratic linear trend gof(fit.jt.quad) rbind(gof(fit.jt),gof(fit.jt.bs),gof(fit.jt.linear),gof(fit.jt.quad)) # Estimates, standard error and Z-statistics require(xtable) temp1 <- summary(fit.jt, print = "Regression", std.error = TRUE, response = 1) round(temp1[[1]][[1]],3) temp2 <- summary(fit.jt, print = "Regression", std.error = TRUE, response = 2) round(temp2[[2]][[1]],3) # Fitted values -------------------------------------------------------- Alt <- as.factor(1:5) Sex = c("F", "M") Method = c("Trampa","Escopeta") Month = 1:33 newdata <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata$Alt) <- levels(data$Alt) levels(newdata$Sex) <- levels(data$Sex) levels(newdata$Method) <- levels(data$Method) newdata$resp <- 1 newdata2 <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata2$Sex) <- levels(data$Sex) levels(newdata2$Method) <- levels(data$Method) levels(newdata$Alt) <- levels(data$Alt) newdata2$resp <- 2 X_new1 <- model.matrix(~ Method*Alt + Sex + Alt*poly(Month, 3), data=newdata) X_new2 <- model.matrix(~ Method*Alt + Sex*Alt + Alt*poly(Month, 4), data = newdata2) names(newdata2) <- names(newdata) fullnewdata <- rbind(newdata, newdata2) levels(fullnewdata$Sex) <- c("Female","Male") levels(fullnewdata$Alt) <- c("Alt-1","Alt-2","Alt-3","Alt-4","Alt-5") levels(fullnewdata$Method) <- c("Firearm","Snare") X <- bdiag(X_new1, X_new2) pred <- mc_link_function(beta = coef(fit.jt, type = "beta")$Estimates, X = X, offset = log(1), link = "log") pontual <- X%*%coef(fit.jt, type = "beta")$Estimates nbeta <- length(coef(fit.jt, type = "beta")$Estimates) VCOV = X%*%vcov(fit.jt)[1:nbeta,1:nbeta]%*%t(X) std.error <- sqrt(diag(VCOV)) pred.min <- pontual - qnorm(0.975)*std.error pred.max <- pontual + qnorm(0.975)*std.error fullnewdata$pred <- pred$mu fullnewdata$pred.min <- exp(pred.min) fullnewdata$pred.max <- exp(pred.max) fullnewdata$FACT <- paste(fullnewdata$Sex,fullnewdata$Method, sep = ":") data$FACT <- paste(data$Sex,data$Method) source("PlotFunctions.R") p.bd <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 10), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 1), type = "l")) + as.layer(xyplot(BD/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.bd) p.ot <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 5), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 2), type = "l")) + as.layer(xyplot(OT/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.ot) ##### Splines model # Fitted values -------------------------------------------------------- Alt <- as.factor(1:5) Sex = c("F", "M") Method = c("Trampa","Escopeta") Month = 1:33 newdata <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata$Alt) <- levels(data$Alt) levels(newdata$Sex) <- levels(data$Sex) levels(newdata$Method) <- levels(data$Method) newdata$resp <- 1 newdata2 <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata2$Sex) <- levels(data$Sex) levels(newdata2$Method) <- levels(data$Method) levels(newdata$Alt) <- levels(data$Alt) newdata2$resp <- 2 X_new1 <- model.matrix(~ Method*Alt + Sex + Alt*bs(Month, 3), data=newdata) X_new2 <- model.matrix(~ Method*Alt + Sex*Alt + Alt*bs(Month, 4), data = newdata2) names(newdata2) <- names(newdata) fullnewdata <- rbind(newdata, newdata2) levels(fullnewdata$Sex) <- c("Female","Male") levels(fullnewdata$Alt) <- c("Alt-1","Alt-2","Alt-3","Alt-4","Alt-5") levels(fullnewdata$Method) <- c("Firearm","Snare") X <- bdiag(X_new1, X_new2) pred <- mc_link_function(beta = coef(fit.jt.bs, type = "beta")$Estimates, X = X, offset = log(1), link = "log") pontual <- X%*%coef(fit.jt.bs, type = "beta")$Estimates nbeta <- length(coef(fit.jt.bs, type = "beta")$Estimates) VCOV = X%*%vcov(fit.jt.bs)[1:nbeta,1:nbeta]%*%t(X) std.error <- sqrt(diag(VCOV)) pred.min <- pontual - qnorm(0.975)*std.error pred.max <- pontual + qnorm(0.975)*std.error fullnewdata$pred <- pred$mu fullnewdata$pred.min <- exp(pred.min) fullnewdata$pred.max <- exp(pred.max) fullnewdata$FACT <- paste(fullnewdata$Sex,fullnewdata$Method, sep = ":") data$FACT <- paste(data$Sex,data$Method) source("PlotFunctions.R") p.bd <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 10), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 1), type = "l")) + as.layer(xyplot(BD/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.bd) p.ot <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 5), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 2), type = "l")) + as.layer(xyplot(OT/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.ot) ### LINEAR TREND # Fitted values -------------------------------------------------------- Alt <- as.factor(1:5) Sex = c("F", "M") Method = c("Trampa","Escopeta") Month = 1:33 newdata <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata$Alt) <- levels(data$Alt) levels(newdata$Sex) <- levels(data$Sex) levels(newdata$Method) <- levels(data$Method) newdata$resp <- 1 newdata2 <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata2$Sex) <- levels(data$Sex) levels(newdata2$Method) <- levels(data$Method) levels(newdata$Alt) <- levels(data$Alt) newdata2$resp <- 2 X_new1 <- model.matrix(~ Method*Alt + Sex + Alt*Month, data=newdata) X_new2 <- model.matrix(~ Method*Alt + Sex*Alt + Alt*Month, data = newdata2) names(newdata2) <- names(newdata) fullnewdata <- rbind(newdata, newdata2) levels(fullnewdata$Sex) <- c("Female","Male") levels(fullnewdata$Alt) <- c("Alt-1","Alt-2","Alt-3","Alt-4","Alt-5") levels(fullnewdata$Method) <- c("Firearm","Snare") X <- bdiag(X_new1, X_new2) pred <- mc_link_function(beta = coef(fit.jt.linear, type = "beta")$Estimates, X = X, offset = log(1), link = "log") pontual <- X%*%coef(fit.jt.linear, type = "beta")$Estimates nbeta <- length(coef(fit.jt.linear, type = "beta")$Estimates) VCOV = X%*%vcov(fit.jt.linear)[1:nbeta,1:nbeta]%*%t(X) std.error <- sqrt(diag(VCOV)) pred.min <- pontual - qnorm(0.975)*std.error pred.max <- pontual + qnorm(0.975)*std.error fullnewdata$pred <- pred$mu fullnewdata$pred.min <- exp(pred.min) fullnewdata$pred.max <- exp(pred.max) fullnewdata$FACT <- paste(fullnewdata$Sex,fullnewdata$Method, sep = ":") data$FACT <- paste(data$Sex,data$Method) source("PlotFunctions.R") p.bd <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 10), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 1), type = "l")) + as.layer(xyplot(BD/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.bd) p.ot <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 5), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 2), type = "l")) + as.layer(xyplot(OT/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.ot) # Quadratic trend ------------------------------------------------------ Alt <- as.factor(1:5) Sex = c("F", "M") Method = c("Trampa","Escopeta") Month = 1:33 newdata <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata$Alt) <- levels(data$Alt) levels(newdata$Sex) <- levels(data$Sex) levels(newdata$Method) <- levels(data$Method) newdata$resp <- 1 newdata2 <- expand.grid("Sex" = Sex, "Method" = Method, "Alt" = Alt, "Month" = Month) levels(newdata2$Sex) <- levels(data$Sex) levels(newdata2$Method) <- levels(data$Method) levels(newdata$Alt) <- levels(data$Alt) newdata2$resp <- 2 X_new1 <- model.matrix(~ Method*Alt + Sex + Alt*poly(Month,2), data=newdata) X_new2 <- model.matrix(~ Method*Alt + Sex*Alt + Alt*poly(Month,4), data = newdata2) names(newdata2) <- names(newdata) fullnewdata <- rbind(newdata, newdata2) levels(fullnewdata$Sex) <- c("Female","Male") levels(fullnewdata$Alt) <- c("Alt-1","Alt-2","Alt-3","Alt-4","Alt-5") levels(fullnewdata$Method) <- c("Firearm","Snare") X <- bdiag(X_new1, X_new2) pred <- mc_link_function(beta = coef(fit.jt.quad, type = "beta")$Estimates, X = X, offset = log(1), link = "log") pontual <- X%*%coef(fit.jt.quad, type = "beta")$Estimates nbeta <- length(coef(fit.jt.quad, type = "beta")$Estimates) VCOV = X%*%vcov(fit.jt.quad)[1:nbeta,1:nbeta]%*%t(X) std.error <- sqrt(diag(VCOV)) pred.min <- pontual - qnorm(0.975)*std.error pred.max <- pontual + qnorm(0.975)*std.error fullnewdata$pred <- pred$mu fullnewdata$pred.min <- exp(pred.min) fullnewdata$pred.max <- exp(pred.max) fullnewdata$FACT <- paste(fullnewdata$Sex,fullnewdata$Method, sep = ":") data$FACT <- paste(data$Sex,data$Method) source("PlotFunctions.R") p.bd <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 10), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 1), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 1), type = "l")) + as.layer(xyplot(BD/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.bd) p.ot <- xyplot(pred ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), ylab = "Fitted values", lty = 1 , auto.key = TRUE, type = "l", ylim = c(-0.5, 5), layout = c(5,4), par.settings = ps) + as.layer(xyplot(pred.min ~ Month | Alt + FACT, data = subset(fullnewdata, resp == 2), lty = 2, type = "l")) + as.layer(xyplot(pred.max ~ Month | Alt + FACT, lty = 2, data = subset(fullnewdata, resp == 2), type = "l")) + as.layer(xyplot(OT/Offset ~ Month | Alt + FACT, data = data, cex = 0.4, pch = 19)) useOuterStrips(p.ot) # Generic function to compute correlation ------------------------------ compute.correlation <- function(tau1, tau2, VCOV) { rho <- tau2/(tau1 + tau2) D.rho.tau1 <- -tau2/(tau1 + tau2)^2 D.rho.tau2 <- 1/(tau1 + tau2) + D.rho.tau1 D.rho <- c(D.rho.tau1, D.rho.tau2) rho - qnorm(0.975)*sqrt(t(D.rho)%*%VCOV%*%D.rho) rho + qnorm(0.975)*sqrt(t(D.rho)%*%VCOV%*%D.rho) return(c(as.numeric(rho), as.numeric(sqrt(t(D.rho)%*%VCOV%*%D.rho)))) } # Computing the correlation induced by each component # of the matrix linear predictor tau.bd <- coef(fit.jt, type = "tau", response = 1)$Estimates tau1 <- tau.bd[1] tau2 <- tau.bd[2] VCOV.HM = vcov(fit.jt)[c("tau11","tau12"),c("tau11","tau12")] compute.correlation(tau1 = tau.bd[1], tau2 = tau.bd[2], VCOV = VCOV.HM) VCOV.ME = vcov(fit.jt)[c("tau11","tau13"),c("tau11","tau13")] compute.correlation(tau1 = tau.bd[1], tau2 = tau.bd[3], VCOV = VCOV.ME) VCOV.LONG = vcov(fit.jt)[c("tau11","tau14"),c("tau11","tau14")] compute.correlation(tau1 = tau.bd[1], tau2 = tau.bd[4], VCOV = VCOV.LONG) tau.ot <- coef(fit.jt, type = "tau", response = 2)$Estimates VCOV.HM.ot = vcov(fit.jt)[c("tau21","tau22"),c("tau21","tau22")] compute.correlation(tau1 = tau.ot[1], tau2 = tau.ot[2], VCOV = VCOV.HM.ot) save.image("FinalResult12082016.RData") #-----------------------------------------------------------------------