Nothing
plot_AtCtEt <- function(AtCtEt, boot=FALSE, xlab, ylab, main, col, legend)
{
# if(class(AtCtEt)!='AtCtEt_model')
if(!is(AtCtEt,'AtCtEt_model'))
{
stop('The first parameter must be an object obtained from the AtCtEt function.')
}
if(boot==TRUE)
{
if(is.null(AtCtEt$boot)==TRUE)
{
stop('Please first run the AtCtEt model with bootstrapping.')
}
}
model_cur <- AtCtEt
l_a <- model_cur$n_beta_a
l_c <- model_cur$n_beta_c
l_e <- model_cur$n_beta_e
#pheno_m <- c(t(data_m[,1:2]))
#pheno_d <- c(t(data_d[,1:2]))
#T_m <- rep(data_m[,3], each=2)
#T_d <- rep(data_d[,3], each=2)
n_beta <- 1:(l_a+l_c+l_e)
order <- 3
x <- seq(from=model_cur$min_t, to=model_cur$max_t, length.out=500)
if(model_cur$n_beta_a>1)
{
bb_a <- splineDesign(model_cur$knots_a, x = x, ord=order, outer.ok = TRUE)
}else{
bb_a <- splineDesign(model_cur$knots_a, x = x, ord=1, outer.ok = TRUE)
if(model_cur$beta_a[1]==-Inf)
{n_beta[1]=0}
}
if(model_cur$n_beta_c>1)
{
bb_c <- splineDesign(model_cur$knots_c, x = x, ord=order, outer.ok = TRUE)
}else{
bb_c <- splineDesign(model_cur$knots_c, x = x, ord=1, outer.ok = TRUE)
if(model_cur$beta_c[1]==-Inf)
{n_beta[model_cur$n_beta_a+1]=0}
}
if(model_cur$n_beta_e>1)
{
bb_e <- splineDesign(model_cur$knots_e, x = x, ord=order, outer.ok = TRUE)
}else{
bb_e <- splineDesign(model_cur$knots_e, x = x, ord=1, outer.ok = TRUE)
}
points_a <- exp(bb_a%*%model_cur$beta_a)
points_c <- exp(bb_c%*%model_cur$beta_c)
points_e <- exp(bb_e%*%model_cur$beta_e)
#fisher <- solve(model_cur$hessian[2:(1+l_a+l_c),2:(1+l_a+l_c)])
max_v <- max(points_c, points_a, points_e)*1.2
plot(range(x), c(0,max_v), type = "n", xlab = xlab, ylab = ylab, main = main)
index <- 1
# bb <- splineDesign(model_cur$knots_a, x = x, ord=order, outer.ok = TRUE)
lines(x, points_a, col = col[1], lwd = 2)
lines(x, points_c, col = col[2], lwd = 2)
lines(x, points_e, col = col[3], lwd = 2)
if(boot == TRUE)
{
lines(AtCtEt$boot$x, AtCtEt$boot$lower.ci_a, col = "orange" ,lty = 2 , lwd = 0.6)
lines(AtCtEt$boot$x, AtCtEt$boot$upper.ci_a, col = "orange" ,lty = 2 , lwd = 0.6)
polygon(c(AtCtEt$boot$x, rev(AtCtEt$boot$x)),c(AtCtEt$boot$upper.ci_a, rev(AtCtEt$boot$lower.ci_a)),col='grey',border = NA, lty=3, density=20)
lines(AtCtEt$boot$x, AtCtEt$boot$lower.ci_c, col = "green" ,lty = 2 , lwd = 0.6)
lines(AtCtEt$boot$x, AtCtEt$boot$upper.ci_c, col = "green" ,lty = 2 , lwd = 0.6)
polygon(c(AtCtEt$boot$x, rev(AtCtEt$boot$x)),c(AtCtEt$boot$upper.ci_c, rev(AtCtEt$boot$lower.ci_c)),col='grey',border = NA, lty=3, density=20)
lines(AtCtEt$boot$x, AtCtEt$boot$lower.ci_e, col = "yellow" ,lty = 2 , lwd = 0.6)
lines(AtCtEt$boot$x, AtCtEt$boot$upper.ci_e, col = "yellow" ,lty = 2 , lwd = 0.6)
polygon(c(AtCtEt$boot$x, rev(AtCtEt$boot$x)),c(AtCtEt$boot$upper.ci_e, rev(AtCtEt$boot$lower.ci_e)),col='grey',border = NA, lty=3, density=20)
}else{
fisher <- solve(model_cur$hessian[n_beta,n_beta])
if((l_a>1)|(model_cur$beta_a[1]!=-Inf))
{
fisher_a <- fisher[index:l_a,index:l_a]
lower <- rep(NA, length(x))
upper <- rep(NA, length(x))
sd <- rep(NA, length(x))
flag <- 0
for(i in 1:length(x))
{
delta <- t(bb_a[i,])%*%fisher_a%*%bb_a[i,]
if(delta>=0)
{
sd[i] <- sqrt(delta)
lower[i] <- sum(bb_a[i,]*model_cur$beta_a) - 1.96*sd[i]
upper[i] <- sum(bb_a[i,]*model_cur$beta_a) + 1.96*sd[i]
}else{flag <- 1}
}
#if(l_a>1)
#{
lower <- exp(lower)
upper <- exp(upper)
#}
lower <- ifelse(lower<0, 0, lower)
upper <- ifelse(upper>max_v, max_v, upper)
index <- index + l_a
if(flag == 0)
{
lines(x, lower, col = "orange" ,lty = 2 , lwd = 0.6)
lines(x, upper, col = "orange" ,lty = 2 , lwd = 0.6)
polygon(c(x, rev(x)),c(upper, rev(lower)),col='grey',border = NA, lty=3, density=20)
}else{
warning('The variance of one of the estimates from the Delta method for the A component is negative. Please try the bootstrap method for the confidence interval or use a different model.')
}
}
if((l_c>1)|(model_cur$beta_c[1]!=-Inf))
{
fisher_c <- fisher[index:(index+l_c-1),index:(index+l_c-1)]
lower <- rep(NA, length(x))
upper <- rep(NA, length(x))
sd <- rep(NA, length(x))
flag <- 0
for(i in 1:length(x))
{
delta <- t(bb_c[i,])%*%fisher_c%*%bb_c[i,]
if(delta>=0)
{
sd[i] <- sqrt(delta)
lower[i] <- sum(bb_c[i,]*model_cur$beta_c) - 1.96*sd[i]
upper[i] <- sum(bb_c[i,]*model_cur$beta_c) + 1.96*sd[i]
}else{flag <- 1}
}
#if(l_c>1)
#{
lower <- exp(lower)
upper <- exp(upper)
#}
lower <- ifelse(lower<0, 0, lower)
upper <- ifelse(upper>max_v, max_v, upper)
index <- index + l_c
if(flag == 0)
{
lines(x, lower, col = "green" ,lty = 2 , lwd = 0.6)
lines(x, upper, col = "green" ,lty = 2 , lwd = 0.6)
polygon(c(x, rev(x)),c(upper, rev(lower)),col='grey',border = NA, lty=3, density=20)
}else{warning('The variance of one of the estimates from the Delta method for the C component is negative. Please try the bootstrap method for the confidence interval or use a different model.')}
}
fisher_e <- fisher[index:(index+l_e-1),index:(index+l_e-1)]
lower <- rep(NA, length(x))
upper <- rep(NA, length(x))
sd <- rep(NA, length(x))
for(i in 1:length(x))
{
sd[i] <- sqrt(t(bb_e[i,])%*%fisher_e%*%bb_e[i,])
lower[i] <- sum(bb_e[i,]*model_cur$beta_e) - 1.96*sd[i]
upper[i] <- sum(bb_e[i,]*model_cur$beta_e) + 1.96*sd[i]
}
#if(l_e>1)
#{
lower <- exp(lower)
upper <- exp(upper)
#}
lower <- ifelse(lower<0, 0, lower)
upper <- ifelse(upper>max_v, max_v, upper)
lines(x, lower, col = "yellow" ,lty = 2 , lwd = 0.6)
lines(x, upper, col = "yellow" ,lty = 2 , lwd = 0.6)
polygon(c(x, rev(x)),c(upper, rev(lower)),col='grey',border = NA, lty=3, density=20)
} # boot == FALSE
if(legend==TRUE)
{
legend(x[1], max_v, c('Additive genetic component','Common environmental component', 'Unique environmental component'), col = col, lty=c(1,1,1), lwd=c(2,2,2))
}
}
plot_AtDtEt <- function(AtCtEt, boot=FALSE, xlab, ylab, main, col, legend)
{
# if(class(AtCtEt)!='AtCtEt_model')
if(!is(AtCtEt,'AtDtEt_model'))
{
stop('The first parameter must be an object obtained from the AtDtEt function.')
}
if(boot==TRUE)
{
if(is.null(AtCtEt$boot)==TRUE)
{
stop('Please first run the AtDtEt model with bootstrapping.')
}
}
model_cur <- AtCtEt
l_a <- model_cur$n_beta_a
l_d <- model_cur$n_beta_d
l_e <- model_cur$n_beta_e
n_beta <- 1:(l_a+l_d+l_e)
order <- 3
x <- seq(from=model_cur$min_t, to=model_cur$max_t, length.out=500)
if(model_cur$n_beta_a>1)
{
bb_a <- splineDesign(model_cur$knots_a, x = x, ord=order, outer.ok = TRUE)
}else{
bb_a <- splineDesign(model_cur$knots_a, x = x, ord=1, outer.ok = TRUE)
if(model_cur$beta_a[1]==-Inf)
{n_beta[1]=0}
}
if(model_cur$n_beta_d>1)
{
bb_d <- splineDesign(model_cur$knots_d, x = x, ord=order, outer.ok = TRUE)
}else{
bb_d <- splineDesign(model_cur$knots_d, x = x, ord=1, outer.ok = TRUE)
if(model_cur$beta_d[1]==-Inf)
{n_beta[model_cur$n_beta_a+1]=0}
}
if(model_cur$n_beta_e>1)
{
bb_e <- splineDesign(model_cur$knots_e, x = x, ord=order, outer.ok = TRUE)
}else{
bb_e <- splineDesign(model_cur$knots_e, x = x, ord=1, outer.ok = TRUE)
}
points_a <- exp(bb_a%*%model_cur$beta_a)
points_d <- exp(bb_d%*%model_cur$beta_d)
points_e <- exp(bb_e%*%model_cur$beta_e)
#fisher <- solve(model_cur$hessian[2:(1+l_a+l_c),2:(1+l_a+l_c)])
max_v <- max(points_d, points_a, points_e)*1.2
plot(range(x), c(0,max_v), type = "n", xlab = xlab, ylab = ylab, main = main)
index <- 1
# bb <- splineDesign(model_cur$knots_a, x = x, ord=order, outer.ok = TRUE)
lines(x, points_a, col = col[1], lwd = 2)
lines(x, points_d, col = col[2], lwd = 2)
lines(x, points_e, col = col[3], lwd = 2)
if(boot == TRUE)
{
lines(AtCtEt$boot$x, AtCtEt$boot$lower.ci_a, col = "orange" ,lty = 2 , lwd = 0.6)
lines(AtCtEt$boot$x, AtCtEt$boot$upper.ci_a, col = "orange" ,lty = 2 , lwd = 0.6)
polygon(c(AtCtEt$boot$x, rev(AtCtEt$boot$x)),c(AtCtEt$boot$upper.ci_a, rev(AtCtEt$boot$lower.ci_a)),col='grey',border = NA, lty=3, density=20)
lines(AtCtEt$boot$x, AtCtEt$boot$lower.ci_d, col = "green" ,lty = 2 , lwd = 0.6)
lines(AtCtEt$boot$x, AtCtEt$boot$upper.ci_d, col = "green" ,lty = 2 , lwd = 0.6)
polygon(c(AtCtEt$boot$x, rev(AtCtEt$boot$x)),c(AtCtEt$boot$upper.ci_d, rev(AtCtEt$boot$lower.ci_d)),col='grey',border = NA, lty=3, density=20)
lines(AtCtEt$boot$x, AtCtEt$boot$lower.ci_e, col = "yellow" ,lty = 2 , lwd = 0.6)
lines(AtCtEt$boot$x, AtCtEt$boot$upper.ci_e, col = "yellow" ,lty = 2 , lwd = 0.6)
polygon(c(AtCtEt$boot$x, rev(AtCtEt$boot$x)),c(AtCtEt$boot$upper.ci_e, rev(AtCtEt$boot$lower.ci_e)),col='grey',border = NA, lty=3, density=20)
}else{
fisher <- solve(model_cur$hessian[n_beta,n_beta])
if((l_a>1)|(model_cur$beta_a[1]!=-Inf))
{
fisher_a <- fisher[index:l_a,index:l_a]
lower <- rep(NA, length(x))
upper <- rep(NA, length(x))
sd <- rep(NA, length(x))
flag <- 0
for(i in 1:length(x))
{
delta <- t(bb_a[i,])%*%fisher_a%*%bb_a[i,]
if(delta>=0)
{
sd[i] <- sqrt(delta)
lower[i] <- sum(bb_a[i,]*model_cur$beta_a) - 1.96*sd[i]
upper[i] <- sum(bb_a[i,]*model_cur$beta_a) + 1.96*sd[i]
}else{flag <- 1}
}
#if(l_a>1)
#{
lower <- exp(lower)
upper <- exp(upper)
#}
lower <- ifelse(lower<0, 0, lower)
upper <- ifelse(upper>max_v, max_v, upper)
index <- index + l_a
if(flag == 0)
{
lines(x, lower, col = "orange" ,lty = 2 , lwd = 0.6)
lines(x, upper, col = "orange" ,lty = 2 , lwd = 0.6)
polygon(c(x, rev(x)),c(upper, rev(lower)),col='grey',border = NA, lty=3, density=20)
}else{
warning('The variance of one of the estimates from the Delta method for the A component is negative. Please try the bootstrap method for the confidence interval or use a different model.')
}
}
if((l_d>1)|(model_cur$beta_d[1]!=-Inf))
{
fisher_d <- fisher[index:(index+l_d-1),index:(index+l_d-1)]
lower <- rep(NA, length(x))
upper <- rep(NA, length(x))
sd <- rep(NA, length(x))
flag <- 0
for(i in 1:length(x))
{
delta <- t(bb_d[i,])%*%fisher_d%*%bb_d[i,]
if(delta>=0)
{
sd[i] <- sqrt(delta)
lower[i] <- sum(bb_d[i,]*model_cur$beta_d) - 1.96*sd[i]
upper[i] <- sum(bb_d[i,]*model_cur$beta_d) + 1.96*sd[i]
}else{flag <- 1}
}
#if(l_c>1)
#{
lower <- exp(lower)
upper <- exp(upper)
#}
lower <- ifelse(lower<0, 0, lower)
upper <- ifelse(upper>max_v, max_v, upper)
index <- index + l_d
if(flag == 0)
{
lines(x, lower, col = "green" ,lty = 2 , lwd = 0.6)
lines(x, upper, col = "green" ,lty = 2 , lwd = 0.6)
polygon(c(x, rev(x)),c(upper, rev(lower)),col='grey',border = NA, lty=3, density=20)
}else{warning('The variance of one of the estimates from the Delta method for the D component is negative. Please try the bootstrap method for the confidence interval or use a different model.')}
}
fisher_e <- fisher[index:(index+l_e-1),index:(index+l_e-1)]
lower <- rep(NA, length(x))
upper <- rep(NA, length(x))
sd <- rep(NA, length(x))
for(i in 1:length(x))
{
sd[i] <- sqrt(t(bb_e[i,])%*%fisher_e%*%bb_e[i,])
lower[i] <- sum(bb_e[i,]*model_cur$beta_e) - 1.96*sd[i]
upper[i] <- sum(bb_e[i,]*model_cur$beta_e) + 1.96*sd[i]
}
#if(l_e>1)
#{
lower <- exp(lower)
upper <- exp(upper)
#}
lower <- ifelse(lower<0, 0, lower)
upper <- ifelse(upper>max_v, max_v, upper)
lines(x, lower, col = "yellow" ,lty = 2 , lwd = 0.6)
lines(x, upper, col = "yellow" ,lty = 2 , lwd = 0.6)
polygon(c(x, rev(x)),c(upper, rev(lower)),col='grey',border = NA, lty=3, density=20)
} # boot == FALSE
if(legend==TRUE)
{
legend(x[1], max_v, c('Additive genetic component','Dominance genetic component', 'Unique environmental component'), col = col, lty=c(1,1,1), lwd=c(2,2,2))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.