## packages and data library("simplexreg") library("plotrix") data("sdac", package = "simplexreg") head(sdac, n = 5) data("retinal", package = "simplexreg") ## plotting (Kernel density): rsimplot <- function(n, mu, sigma) { x <- seq(0.01, 0.99, 0.01) for (j in 1:length(sigma)) { for (i in 1:length(mu)) { plot(x, dsimplex(x, mu[i], sigma[j]), ylim = c(0, max(dsimplex(x, mu[i], sigma[j])) + 0.5), type = "l", col = "blue") lines(density(rsimplex(n, mu[i], sigma[j])), lty = 2) } } } par(mfrow = c(3, 3), mar = c(2, 2, 2, 2), cex = 0.55) rsimplot(5000, c(0.1, 0.5, 0.7), c(4, 2, 1)) ## or plotting (histogram): rsimplot <- function(n, mu, sigma) { x <- seq(0.01, 0.99, 0.01) for (j in 1:length(sigma)) { for (i in 1:length(mu)) { hist(rsimplex(n, mu[i], sigma[j]), density = 50, xlim = c(0, 1), ylim = c(0, max(dsimplex(x, mu[i], sigma[j])) + 0.5), freq = F, breaks = 20, xlab = " ", ylab = " ", main = NULL) lines(x, dsimplex(x, mu[i], sigma[j]), type = "l", col = "blue") } } } par(mfrow = c(3, 3), mar = c(2, 2, 2, 2), cex = 0.55) rsimplot(5000, c(0.1, 0.5, 0.7), c(4, 2, 1)) ## PBSC study # homogeneous model sim.glm1 <- simplexreg(rcd ~ ageadj + chemo, data = sdac) summary(sim.glm1) # heterogeneous model sim.glm2 <- simplexreg(rcd ~ ageadj + chemo | age, data = sdac) summary(sim.glm2) AIC(sim.glm1, sim.glm2) # residual analysis plot(sim.glm1, type = "residuals", res = "stdPerr", ylim = c(-3, 3)) plot(sim.glm2, type = "residuals", res = "stdPerr", ylim = c(-3, 3)) ## Eye surgery study # examine the correlation structure data("retinal", package = "simplexreg") sim.glm3 <- simplexreg(Gas ~ LogT + LogT2 + Level | LogT+Level | Time, data = retinal, id = ID) plot(sim.glm3, type = "corr", xlim = c(-2.5, 2.5), ylim = c(-2.5, 2.5), pch = 16) plot(sim.glm3, type = "corr", lag = 2, xlim = c(-2.5, 2.5), ylim = c(-2.5, 2.5), pch = 16) plot(sim.glm3, type = "corr", lag = 3, xlim = c(-2.5, 2.5), ylim = c(-2.5, 2.5), pch = 16) plot(sim.glm3, type = "corr", lag = 4, xlim = c(-2.5, 2.5), ylim = c(-2.5, 2.5), pch = 16) # fit the data sim.gee2 <- simplexreg(Gas ~ LogT + LogT2 + Level | LogT + Level | Time, corr = "AR1", id = ID, data = retinal) summary(sim.gee2) # residual analysis and partial deviance plot(sim.gee2, type = "residuals", ylim = c(-6, 6), pch = 16) plot(sim.gee2, type = "GOF", xlab = "Days after Gas Injection", ylim = c(0, 40))