Telomere Data Analysis
Introduction
blah blah blah…
Preliminaries
# code chunk
# Data for live-coding exercies
# Lowest p-value in the subset of data
#
# (cntrl = control = ambient sounds)
# (dist = disturbed = sound disturbance)
#
# 6 rows of data
# telomere length for control (in kilobases)
# telomere length for sound disturbance (in kilobases)
# corticosterone concentration control
# corticosterone concentration distrubrance
# tretment
# treatment
# How to enter these data into vectors?
# Task 1) The lines of data below need to be formatted for loading into vectors
# name them as follows and load the data into vectors.
# "telos.cntrl"
# "telos.dist"
# "cort.cntrl"
# "cort.dist
# "trt.cntrl"
# "trt.dist
# data in short vector for each trt
telos.cntrl <- c(1.202 ,1.135, 1.116, 1.339 ,0.948 ,1.194 ,1.179)
telos.dist <- c(0.829 ,1.121, 1.128, 1.002, 0.866, 1.087, 0.935, 0.929, 1.012)
cort.cntrl <- c(5.93, 10.45, 1.84, 8.89, 6.32, 1.95, 1.47)
cort.dist <- c(7.38, 0.85, 4.76, 7.26, 5.39, 1.36, 18.35, 5.73, 4.22)
trt.cntrl <- c("cntrl", "cntrl", "cntrl", "cntrl", "cntrl", "cntrl", "cntrl")
trt.dist <- c("dist", "dist", "dist", "dist", "dist", "dist", "dist", "dist", "dist")
# data in long
telos <- c(telos.cntrl, telos.dist)
trt <- c(trt.cntrl, trt.dist)Summary statistics
# means
mean.telos.cntrl <- mean(telos.cntrl)
mean.telos.dist <- mean(telos.dist)
print(mean.telos.cntrl)## [1] 1.159print(mean.telos.dist)## [1] 0.9898889# SD
sd.telos.cntrl <- sd(telos.cntrl)
sd.telos.dist <- sd(telos.dist)
# sample size
n.telos.cntrl <- length(telos.cntrl)
n.telos.dist <- length(telos.dist)
# Standard error (SE)
## SE = SD/sqrt(N)
0.1174876/2.645751## [1] 0.04440614se.telos.cntrl <- sd.telos.cntrl/sqrt(n.telos.cntrl)
se.telos.dist <- sd.telos.dist/sqrt(n.telos.dist)
# Confidence interval
## CI "Arm" = 2*SE
upper.arm.cntrl <- 2*se.telos.cntrl + mean.telos.cntrlThe lowbrow_errbar() function
RUN THIS ENTIRE CODE CHUNK TO LOAD THE FUNCTION MUST RUN THE WHOLE THING _ its big - clikc on the green triangle in the corner of the code chunk
# This is a simple function that plots means as individual points (not bars) and
# approximate 95% confidence intervals using user-supplied means and standard
# errors or 95% CIs. It provides some feed back via the console and the plot
# itself regarding errors and how to make the plot more polished.
#
# Base R has no default function to make plots like this, though there are
# numerous packages that have function that can do this
# (sciplot, ggplot, Hmisc, gplot, psych). I wrote this function as bare-bones,
# "low brow"alternative that my students could load via a script and which has
# only a few, simple arguements. I teach them to calculate means and standard
# deviations using the tapply() function.
#
# I’ve found that loading packages can be difficult on computers behind campus
# firewalls, as can updating R if a package requires an up-to-date installation.
# At the end of this document I provide notes about and links to other ways to
# approach this problem with R published packages.
# For examples and test of the function see
# https://rpubs.com/brouwern/plotmeans2
# means = your means. , contained in a vector, eg c(mean1, means2...)
# SEs = your standard errors, contained in a vector,
# CI.hi = upper confidence intervals, in a vector,
# CI.lo = lower confidence itnerval, in a vector
# categories = the names of the categories/groups, in the order that they appear
# x.axis.label = what should be plotted on the x axis
# y.axis.label = what should be plotted on the y axis
# adjust.y.max = allows you to adjust the y axis
# adjust.y.min
# adjust.x.spacing
#### The function STARTS here ####
lowbrow_errbars <- function(means = NULL,
SEs = NULL,
CI.hi = NULL,
CI.lo = NULL,
categories = NULL,
x.axis.label = "Groups",
y.axis.label = "'y.axis.label' sets the axis label",
adjust.y.max = 0,
adjust.y.min = 0,
adjust.x.spacing = 5){
# Error messages
if(is.null(means)==T){
stop("No means entered") }
if(is.null(SEs)==T & is.null(CI.lo) == T){
stop("No standard errors entered") }
if(is.null(CI.lo) != T & is.null(CI.hi)){
stop("CI.lo entered but no CI.hi") }
if(is.null(CI.hi) != T & is.null(CI.lo)){
stop("CI.hi entered but no CI.lo") }
#check if both SE and and CI.is enter
if(is.null(SEs) == F & is.null(CI.hi) == F){
stop("Both SEs and CI.hi entered; use eithe SEs or CIs") }
if(is.null(SEs) == F & is.null(CI.lo) == F){
stop("Both SEs and CI.lo entered; use eithe SEs or CIs") }
#Check if the number of means matches the number of SE
n.means <- length(means)
n.SEs <- length(SEs)
n.CI.lo <- length(CI.lo)
n.CI.hi <- length(CI.hi)
#CHeck standard errors against means
if(n.means != n.SEs & is.null(CI.lo)==T){
error.message <- paste("The number of means is",n.means,"but the number of standard errors is",n.SEs,sep = " ")
stop(error.message) }
#Check CIs against CIs
if(n.CI.lo != n.CI.hi){
error.message <- paste("The number of CI.lo is",n.CI.lo,"but the number of standard errors is",n.CI.hi,sep = " ")
stop(error.message) }
#assign arbitrary categories
if(is.null(categories) == T) {
categories <- paste("Group",1:n.means)
print("Set categoris labls with 'categories=' ")
}
# calculate values for plotting limits
if(is.null(SEs)==F){
y.max <- max(means+2*SEs) + adjust.y.max
y.min <- min(means-2*SEs) - adjust.y.min
}
# calculate values for plotting limits
if(is.null(SEs)==T){
y.max <- max(CI.hi) + adjust.y.max
y.min <- min(CI.lo) - adjust.y.min
}
#determine where to plot points along x-axis
x.values <- 1:n.means
x.values <- x.values/adjust.x.spacing
#set x axis min/max
x.axis.min <- min(x.values)-0.05
x.axis.max <- max(x.values)+0.05
x.limits <- c(x.axis.min,x.axis.max)
#Plot means
plot(means ~ x.values,
xlim = x.limits,
ylim = c(y.min,y.max),
xaxt = "n",
xlab = "",
ylab = "",
cex = 1.25,
pch = 16)
#Add x labels
axis(side = 1,
at = x.values,
labels = categories
)
# Plot confidence intervals
if(is.null(CI.hi) == FALSE &
is.null(CI.lo) == FALSE){
#Plot upper error bar for CIs
lwd. <- 2
arrows(y0 = means,
x0 = x.values,
y1 = CI.hi,
x1 = x.values,
length = 0,
lwd = lwd.)
#Plot lower error bar
arrows(y0 = means,
x0 = x.values,
y1 = CI.lo,
x1 = x.values,
length = 0,
lwd = lwd.)
}
# Estimate CIs from SEs
if(is.null(SEs) == FALSE &
is.null(CI.lo) == TRUE){
lwd. <- 2
arrows(y0 = means,
x0 = x.values,
y1 = means+2*SEs,
x1 = x.values,
length = 0,
lwd = lwd.)
#Plot lower error bar
arrows(y0 = means,
x0 = x.values,
y1 = means-2*SEs,
x1 = x.values,
length = 0,
lwd = lwd.)
}
mtext(text = x.axis.label,side = 1,line = 2)
mtext(text = y.axis.label,side = 2,line = 2)
mtext(text = "Error bars = 95% CI",side = 3,line = 0,adj = 0)
}
#### The function ENDS here ####
#### The function ENDS here ####
#### The function ENDS here ####
#### The function ENDS here ####
#### The function ENDS here ####
#
# Other approaches
# These approach vary in if/how they generate the intial plot of the means before adding the error bars.
#
# gplot::plotmeans http://svitsrv25.epfl.ch/R-doc/library/gplots/html/plotmeans.html
#
# Using ggplot w/geom_errorbar ggplot is the standard now for advanced graphics. http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
#
# https://www.r-bloggers.com/line-plot-for-two-way-designs-using-ggplot2/
#
# Using sciplot:lineplot.CI See examples above
#
# Hmisc::errbar http://svitsrv25.epfl.ch/R-doc/library/Hmisc/html/errbar.html
#
# Using psych::error.bars https://www.personality-project.org/r/html/error.bars.html
# If you’re into “dynamite plots” https://www.r-bloggers.com/building-barplots-with-error-bars/ http://monkeysuncle.stanford.edu/?p=485Use the lowbrow_errbar() function
# don't ==
# don't <-
lowbrow_errbars(means = c(mean.telos.cntrl, # means for 2 treatments
mean.telos.dist),
SEs = c(se.telos.cntrl, # SEs for 2 treatments
se.telos.dist),
categories = c("Control","Disturbance"),
x.axis.label = "Experimental treatment",
y.axis.label = "Telomere length (kb)"
)
#make a dataframe
df <- data.frame(telos, trt)# t test
# y = x
# y = f(x)
# y ~ x
# telo ~ treatment
t.test(telos ~ trt, data = df)##
## Welch Two Sample t-test
##
## data: telos by trt
## t = 2.9522, df = 12.485, p-value = 0.01165
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.04483683 0.29338539
## sample estimates:
## mean in group cntrl mean in group dist
## 1.1590000 0.9898889