library("BAMBI") # ****************************************************************** # plot wnorm density # ****************************************************************** mu <- pi kappa <- c(0.01, 1, 10) phi <- seq(0, 2*pi, length.out = 100) dens <- sapply(kappa, function(j) dwnorm(x = phi, mu = mu, kappa = j)) cols <- c("orange", "blue", "grey") pdf("wnormden_model.pdf") par(mar = c(5, 5, 1, 1)) plot(phi, dens[, 1], type = "l", ylab = bquote("f"["WN"]*"("*phi~"|"~mu*", "*kappa*")"), xlab = expression(phi), ylim = range(dens), col = cols[1], cex.axis = 1.5, cex.lab = 1.5, xaxt = "n") axis(side = 1, at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi)), cex.axis = 1.5, cex.lab = 2) for (i in 2:3) { points(phi, dens[, i], col = cols[i], type = "l") } legend_expr <- parse(text = paste0("mu*\' = \'*pi*\', \'*", "kappa*\' = \'*", kappa)) legend("topright", legend = legend_expr, col = cols, lty = 1, cex = 1.3) dev.off() # ****************************************************************** # plot vm density # ****************************************************************** mu <- pi kappa <- c(0.01, 1, 10) phi <- seq(0, 2*pi, length.out = 100) dens <- sapply(kappa, function(j) dvm(x = phi, mu = mu, kappa = j)) cols <- c("orange", "blue", "grey") pdf("vmden_model.pdf") par(mar = c(5, 5, 1, 1)) plot(phi, dens[, 1], type = "l", ylab = bquote("f"["vM"]*"("*phi~"|"~mu*", "*kappa*")"), xlab = expression(phi), ylim = range(dens), col = cols[1], xaxt = "n", cex.axis = 1.5, cex.lab = 1.5) axis(side = 1, at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi)), cex.axis = 1.5, cex.lab = 2) for (i in 2:3) { points(phi, dens[, i], col = cols[i], type = "l") } legend_expr <- parse(text = paste0("mu*\' = \'*pi*\', \'*", "kappa*\' = \'*", kappa)) legend("topright", legend = legend_expr, col = cols, lty = 1, cex = 1.3) dev.off() # ****************************************************************** # plot wnorm2 density # ****************************************************************** mu1 <- mu2 <- pi kappa1 <- kappa2 <- c(0.01, 0.5, 2, 1, 1, 1) kappa3 <- c(0, 0, 0, -0.5, 0, 0.5) phi <- psi <- seq(0, 2*pi, length.out = 50) dens <- sapply(seq(length(kappa3)), function(j) dwnorm2(x = cbind(phi, psi), mu1 = mu1, mu2 = mu2, kappa1 = kappa1[j], kappa2 = kappa2[j], kappa3 = kappa3[j])) dens_surf <- lapply(seq(length(kappa3)), function(j) surface_model("wnorm2", mu1 = mu1, mu2 = mu2, kappa1 = kappa1[j], kappa2 = kappa2[j], kappa3 = kappa3[j], xpoints = phi, ypoints = psi, xlab = expression(phi), ylab = expression(psi), zlim = range(dens), # main = paste0("mu1 = mu2 = ", expression(pi), # "\nkappa1 = kappa2 = ", kappa1[j], # "\nkappa3 = ", kappa3[j]), main = bquote(kappa[1] == .(kappa1[j])* "," ~ kappa[2] == .(kappa1[j])* "," ~ kappa[3] == .(kappa3[j])), scales = list( arrows = FALSE, x = list(at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi))), y = list(at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi))), col = "black", font = 1, tck = 1))) pdf("wnorm2den_model.pdf", height = 13, width = 19) print(gridExtra::grid.arrange(grobs = dens_surf, ncol = 3)) dev.off() # ****************************************************************** # plot vmsin density # ****************************************************************** mu1 <- mu2 <- pi kappa1 <- kappa2 <- 1 kappa3 <- c(-2, -0.5, 0, 0.5, 1, 2) phi <- psi <- seq(0, 2*pi, length.out = 50) dens <- sapply(seq(length(kappa3)), function(j) dvmsin(x = cbind(phi, psi), mu1 = mu1, mu2 = mu2, kappa1 = kappa1, kappa2 = kappa2, kappa3 = kappa3[j])) dens_surf <- lapply(seq(length(kappa3)), function(j) surface_model("vmsin", mu1 = mu1, mu2 = mu2, kappa1 = kappa1, kappa2 = kappa2, kappa3 = kappa3[j], xpoints = phi, ypoints = psi, xlab = expression(phi), ylab = expression(psi), zlim = range(dens), # main = paste0("mu1 = mu2 = ", expression(pi), # "\nkappa1 = kappa2 = ", kappa1, # "\nkappa3 = ", kappa3[j]), main = bquote(kappa[1] == .(kappa1)* "," ~ kappa[2] == .(kappa1)* "," ~ kappa[3] == .(kappa3[j])), scales = list( arrows = FALSE, x = list(at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi))), y = list(at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi))), col = "black", font = 1, tck=1))) pdf("vmsinden_model.pdf", height = 13, width = 19) print(gridExtra::grid.arrange(grobs = dens_surf, ncol = 3)) dev.off() # ****************************************************************** # plot vmcos density # ****************************************************************** mu1 <- mu2 <- pi kappa1 <- kappa2 <- 1 kappa3 <- c(-2, -0.5, 0, 0.5, 1, 2) phi <- psi <- seq(0, 2*pi, length.out = 50) dens <- sapply(seq(length(kappa3)), function(j) dvmcos(x = cbind(phi, psi), mu1 = mu1, mu2 = mu2, kappa1 = kappa1, kappa2 = kappa2, kappa3 = kappa3[j])) dens_surf <- lapply(seq(length(kappa3)), function(j) surface_model("vmcos", mu1 = mu1, mu2 = mu2, kappa1 = kappa1, kappa2 = kappa2, kappa3 = kappa3[j], xpoints = phi, ypoints = psi, xlab = expression(phi), ylab = expression(psi), zlim = range(dens), # main = paste0("mu1 = mu2 = ", expression(pi), # "\nkappa1 = kappa2 = ", kappa1, # "\nkappa3 = ", kappa3[j]), main = bquote(kappa[1] == .(kappa1)* "," ~ kappa[2] == .(kappa1)* "," ~ kappa[3] == .(kappa3[j])), scales = list( arrows = FALSE, x = list(at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi))), y = list(at = c(0, pi/3, 2*pi/3, pi, 4*pi/3, 5*pi/3, 2*pi), label = c("0", expression(pi/3), expression(2*pi/3), expression(pi), expression(4*pi/3), expression(5*pi/3), expression(2*pi))), col = "black", font = 1, tck=1))) pdf("vmcosden_model.pdf", height = 13, width = 19) print(gridExtra::grid.arrange(grobs = dens_surf, ncol = 3)) dev.off() #----------------------------------------------------------------- # Ramachandran (scatter) plot for 8TIM data pdf("rama.pdf") plot(tim8, pch = 16, xlim = c(0, 2*pi), ylim = c(0, 2*pi), main = "8TIM", col = scales::alpha("black", 0.3), xlab = expression(phi), ylab = expression(psi), cex.axis = 1.3, cex.lab = 1.3) dev.off() # ****************************************************************** # Data analysis with vmsin model # ****************************************************************** #----------------------------------------------------------------- # parallelization library("future") plan(multisession, gc = TRUE) #----------------------------------------------------------------- #----------------------------------------------------------------- # fitting 4 comp vmsin set.seed(12321) fit.vmsin.4comp <- fit_angmix("vmsin", tim8, ncomp = 4, n.iter = 2e4, n.chains = 3) #----------------------------------------------------------------- #----------------------------------------------------------------- # fitting optimum mixture of vmsin set.seed(12321) fit.vmsin <- fit_incremental_angmix(model="vmsin", data = tim8, crit = "LOOIC", start_ncomp = 2, max_ncomp = 10, n.iter = 2e4, n.chains = 3) # get the maximum log likelihood of the intermediate fits fit.vmsin$maxllik.all # extract the best model fit.vmsin.best <- bestmodel(fit.vmsin) # lpd and param traces pdf("vmsin_lpdtrace.pdf") lpdtrace(fit.vmsin.best) dev.off() pdf("vmsin_paramtrace.pdf") paramtrace(fit.vmsin.best) dev.off() bitmap("vmsin_lpdtrace.png", width = 7, height = 7, units = 'in', res = 300) lpdtrace(fit.vmsin.best) dev.off() bitmap("vmsin_paramtrace_comp1_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 1) dev.off() bitmap("vmsin_paramtrace_comp2_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 2) dev.off() bitmap("vmsin_paramtrace_comp3_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 3) dev.off() bitmap("vmsin_paramtrace_comp4_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 4) dev.off() # convert to mcmc objects mcmc.vmsin.best <- coda::as.mcmc.list(fit.vmsin.best) geweke_res <- coda::geweke.diag(mcmc.vmsin.best, frac1 = 0.5, frac2 = 0.5) pdf("vmsin_geweke_diag_all.pdf", height = 3.5, width = 6) par(mfrow = c(1, 3), mar = c(5, 6, 2, 1)) for(j in 1:3) { barplot(geweke_res[[j]]$z, horiz = TRUE, names = names(geweke_res[[j]]$z), xlim = range(geweke_res[[j]]$z, -3, 3), ylab = "", xaxt = 'n', las = 2) axis(1, las = 1) title(main = paste("Chain", j), xlab = "geweke.diag") } dev.off() # fix label switchings fit.vmsin.best <- fix_label(fit.vmsin.best) # redraw param traces pdf("vmsin_paramtrace.pdf") paramtrace(fit.vmsin.best) dev.off() bitmap("vmsin_paramtrace_fix_comp1_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 1) dev.off() bitmap("vmsin_paramtrace_fix_comp2_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 2) dev.off() bitmap("vmsin_paramtrace_fix_comp3_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 3) dev.off() bitmap("vmsin_paramtrace_fix_comp4_%d.png", width = 7, height = 7, units = 'in', res = 300) paramtrace(fit.vmsin.best, comp.label = 4) dev.off() # find point estimates round(pointest(fit.vmsin.best, fn = "MODE"), 2) round(pointest(fit.vmsin.best, fn = mean), 2) # contourplot pdf("vmsincontour_map.pdf", height = 5.5, width = 5.5) contour(fit.vmsin.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # density surface pdf("vmsinden_map.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.vmsin.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # contourplot pdf("vmsincontour.pdf", height = 5.5, width = 5.5) contour(fit.vmsin.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # density surface pdf("vmsinden.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.vmsin.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # summary of the estimated parameters summary(fit.vmsin.best) # find likelihood of the fitted model logLik(fit.vmsin.best) # generate from the fitted model set.seed(12321) vmsin.data <- r_fitted(nrow(tim8), fit.vmsin.best, fn = "MODE") pdf("vmsinrama.pdf") plot(vmsin.data, xlab = expression(phi), ylab = expression(psi), cex.axis = 1.3, cex.lab = 1.3, xlim = c(0, 2*pi), pch = 16, ylim = c(0, 2*pi), col = scales::alpha("black", 0.3)) title("Data generated from best fitted vmsin") dev.off() logLik(fit.vmsin.best) set.seed(12321) fit.vmsin$crit.all # ****************************************************************** # vmcos model # ****************************************************************** #----------------------------------------------------------------- # fitting optimum mixture of vmcos set.seed(12321) fit.vmcos <- fit_incremental_angmix(model="vmcos", data = tim8, crit = "LOOIC", start_ncomp = 2, max_ncomp = 10, n.iter = 2e4, n.chains = 3, use_best_chain = FALSE, n_qrnd = 1e4) # get the maximum log likelihood of the intermediate fits fit.vmcos$maxllik.all # extract the best model fit.vmcos.best <- bestmodel(fit.vmcos) # lpd trace lpdtrace(fit.vmcos.best) # parameter trace paramtrace(fit.vmcos.best) # fix label switchings fit.vmcos.best <- fix_label(fit.vmcos.best) # contourplot pdf("vmcoscontour_map.pdf", height = 5.5, width = 5.5) contour(fit.vmcos.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # density surface pdf("vmcosden_map.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.vmcos.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # contourplot pdf("vmcoscontour.pdf", height = 5.5, width = 5.5) contour(fit.vmcos.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # density surface pdf("vmcosden.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.vmcos.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # summary of the estimated parameters summary(fit.vmcos.best) # find likelihood of the fitted model logLik(fit.vmcos.best) fit.vmcos$crit.best # ****************************************************************** # wnorm2 model # ****************************************************************** #----------------------------------------------------------------- # fitting optimum mixture of wnorm2 set.seed(12321) fit.wnorm2 <- fit_incremental_angmix(model="wnorm2", data = tim8, crit = "LOOIC", start_ncomp = 2, max_ncomp = 10, n.iter = 2e4, n.chains = 3, use_best_chain = FALSE) # get the maximum log likelihood of the intermediate fits fit.wnorm2$maxllik.all # extract the best model fit.wnorm2.best <- bestmodel(fit.wnorm2) # lpd and param traces lpdtrace(fit.wnorm2.best) paramtrace(fit.wnorm2.best) # fix label switchings fit.wnorm2.best <- fix_label(fit.wnorm2.best) # contourplot pdf("wnorm2contour_map.pdf", height = 5.5, width = 5.5) contour(fit.wnorm2.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # density surface pdf("wnorm2den_map.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.wnorm2.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # contourplot pdf("wnorm2contour.pdf", height = 5.5, width = 5.5) contour(fit.wnorm2.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # density surface pdf("wnorm2den.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.wnorm2.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # summary of the estimated parameters summary(fit.wnorm2.best) # find likelihood of the fitted model logLik(fit.wnorm2.best) fit.wnorm2$crit.best # ****************************************************************** # compare vmsin, vmcos and wnorm2 fits # ****************************************************************** looic.all.bi <- list( vmsin.4 = bestcriterion(fit.vmsin), vmcos.4 = bestcriterion(fit.vmcos), wnorm2.4 = bestcriterion(fit.wnorm2) ) comp <- loo::loo_compare(x = looic.all.bi) comp find_ci <- function(x, digits = 1) { round(c(lower = unname(x[1] - 2*x[2]), upper = unname(x[1] + 2*x[2])), digits = digits) } t(apply(comp[-1, c("elpd_diff", "se_diff")], 1, find_ci)) #----------------------------------------------------------------- # histogram for wind data pdf("windhist.pdf") hist(wind[, "angle"]) dev.off() # ****************************************************************** # vm model # ****************************************************************** #----------------------------------------------------------------- # fitting optimum mixture of vm #----------------------------------------------------------------- set.seed(12321) fit.vm <- fit_incremental_angmix(model="vm", data = wind[, "angle"], crit = "LOOIC", start_ncomp = 1, max_ncomp = 10, n.iter = 2e4, n.chains = 3) # get the maximum log likelihood of the intermediate fits fit.vm$maxllik.all # extract the best model fit.vm.best <- bestmodel(fit.vm) # # lpd and param traces # pdf("vm_lpdtrace.pdf") # lpdtrace(fit.vm.best) # dev.off() bitmap("vm_lpdtrace.png") lpdtrace(fit.vm.best) dev.off() # pdf("vm_paramtrace.pdf") # paramtrace(fit.vm.best) # dev.off() bitmap("vm_paramtrace_%d.png") paramtrace(fit.vm.best) dev.off() # fix label switchings fit.vm.best <- fix_label(fit.vm.best) round(pointest(fit.vm.best, fn = "MODE"), 2) round(pointest(fit.vm.best, fn = mean), 2) # density curve pdf("vmden_map.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.vm.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # density curve pdf("vmden.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.vm.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # summary of the estimated parameters summary(fit.vm.best) # ****************************************************************** # wnorm model # ****************************************************************** set.seed(12321) fit.wnorm <- fit_incremental_angmix(model = "wnorm", data = wind[, "angle"], crit = "LOOIC", start_ncomp = 1, max_ncomp = 10, n.iter = 2e4, n.chains = 3) # get the maximum log likelihood of the intermediate fits fit.wnorm$maxllik.all # extract the best model fit.wnorm.best <- bestmodel(fit.wnorm) # lpd and param traces bitmap("wnorm_lpdtrace.png") lpdtrace(fit.wnorm.best) dev.off() pdf("wnorm_paramtrace.pdf") paramtrace(fit.wnorm.best) dev.off() # fix label switchings fit.wnorm.best <- fix_label(fit.wnorm.best) # density surface pdf("wnormden_map.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.wnorm.best, fn = "MODE", xlab = expression(phi), ylab = expression(psi)) dev.off() # density curve pdf("wnormden.pdf", height = 5.5, width = 5.5) lattice::densityplot(fit.wnorm.best, fn = mean, xlab = expression(phi), ylab = expression(psi)) dev.off() # summary of the estimated parameters summary(fit.wnorm.best) # ****************************************************************** # compare vm and wnorm model # ****************************************************************** looic.all.uni <- list( vm.2 = bestcriterion(fit.vm), wnorm.2 = bestcriterion(fit.wnorm) ) loo::loo_compare(x = looic.all.uni) sessionInfo()