Telomere Data Analysis
Introductions
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
#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)
telos.dist <- rnorm(n=8, mean = 1000, sd = 30)
telos.cntrl <- rnorm(n=8, mean = 1000, sd = 30)
cort.cntrl5.93 <- c(10.45, 1.84, 8.89, 6.32, 1.95, 1.47)
cort.dist7.38 <- c(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 vectors
telos <- c(telos.cntrl, telos.dist)
trt <- c(trt.cntrl, trt.dist)
# Task 2) How do you calculate the mean for the telomeres? Standard deviation?
# Sample size?Summary Stats
# means
mean.telos.cntrl <- mean(telos.cntrl)
mean.telos.dist <- mean(telos.dist)
# 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)
# std error (SE)
# SE = SD/sqrt(N)
se.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.cntrl
upper.arm.dist <- 2*se.telos.dist + mean.telos.distLoad the R function lowbrow_errbars
#did it get loaded
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)
}
ls()## [1] "cort.cntrl5.93" "cort.dist7.38" "lowbrow_errbars" "mean.telos.cntrl"
## [5] "mean.telos.dist" "n.telos.cntrl" "n.telos.dist" "sd.telos.cntrl"
## [9] "sd.telos.dist" "se.telos.cntrl" "se.telos.dist" "telos"
## [13] "telos.cntrl" "telos.dist" "trt" "trt.cntrl"
## [17] "trt.dist" "upper.arm.cntrl" "upper.arm.dist"exists("lowbrow_errbars")## [1] TRUElowbrow errbars
#dont use == or <- inside functions
lowbrow_errbars(mean = c(mean.telos.cntrl, mean.telos.dist),
SEs = c(se.telos.cntrl, se.telos.dist))## [1] "Set categoris labls with 'categories=' "
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 = 0.3335, df = 10.044, p-value = 0.7456
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -24.29929 32.85985
## sample estimates:
## mean in group cntrl mean in group dist
## 1000.0074 995.7271This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.