Getting help

?FunctionName    #opens a help file on that function
example(x)    #prints an example of the use of function x
R reference card

Reading in Data and viewing

read.delim() #for tab delimited data w/ or w/out header
#example: read.delim("C:/stuff/fickle.csv", header=TRUE)

read.csv()  #this function reads in a comma delimited file with a header
#example :
read.csv("C:/stuff/fickle.csv")

read.table()  #this function reads in a number of different files depending on the arguments
#example:
read.table("C:/stuff/fickle.txt",sep="")
#sep is the argument for the type of delimiting character
#"" is for white space which includes tabs
#\t is for tab
#"," is for comma
#can also specify which columns to coerce to factors

head(x)  #prints top 6 rows of the object x
tail(x)  #prints bottom 6 rows of the object x
str(x)  #prints the structure of the object x

Attaching objects to the search path
attach(data)
detach(data)
search() #to see what's on the search path

Recoding/Recode data example using a simple base statement (definitely most efficient method)
Example 1:
as.factor(dat$date) #make sure working with a factor
temp <- c(37,37,37,37,37,35,33,31,37,35,25,29,27,30,28,29,28,27,27.5,28.5,28.3,27.7)[dat$date]
#temp will now be coded according to temperatures on the given dates
#make sure your vector of temperatures (or whatever) corresponds exactly to the order of the levels
#of the dat$date factor.
#Example 2:
color <- c("green", "blue", "red", "pink")[dat$morph]
#dat$morph is a factor with four levels

Recoding/Recode data example using a for loop and if else statments

#recode morph to lump crescents
dat1 <- data.frame(dat1, morph2 = NA)

for(i in 1:929){
 if (is.na(dat1[i,31])) {dat1[i,50] <- NA}
 else if (dat1[i,31] == "cres-n" |
  dat1[i,31] == "cres-r" |
  dat1[i,31] == "cres-nb"){dat1[i,50] <- "cres"}
 else if (dat1[i,31] == "cush") {dat1[i,50] <- "cush"}
 else if (dat1[i,31] == "flat") {dat1[i,50] <- "flat"}
 else if (dat1[i,31] == "irr") {dat1[i,50] <- "irr"}
 }
dat1$morph2 <- as.factor(dat1$morph2)

Recoding/Recode data example using which() and grep() functions.  More efficient than loops.

#add variables for blocks
#two columns with 5 blocks each (blk_c)
#five rows with 2 blocks each (blk_r)
#individual blocks (10 in total) (blks)

#blk_c
dat <- data.frame(dat, blk_c = 1)
col1 <- which(dat$col == 1) #which returns the row indices that match
col2 <- which(dat$col == 2)
col3 <- which(dat$col == 3)
col4 <- which(dat$col == 4)
col5 <- which(dat$col == 5)
blk_c1 <- c(col1,col2,col3,col4,col5)

col6 <- which(dat$col == 6)
col7 <- which(dat$col == 7)
col8 <- which(dat$col == 8)
col9 <- which(dat$col == 9)
col10 <- which(dat$col == 10)
blk_c2 <- c(col6,col7,col8,col9,col10)

dat$blk_c[blk_c1] <- 1 #here apply the recode value "1" to the matching row indices
dat$blk_c[blk_c2] <- 2

dat$blk_c <- as.factor(dat$blk_c)

#blk_r
dat <- data.frame(dat, blk_r = 1)
row1 <- which(dat$row == 1)
row2 <- which(dat$row == 2)
row3 <- which(dat$row == 3)
row4 <- which(dat$row == 4)
row5 <- which(dat$row == 5)
row6 <- which(dat$row == 6)
row7 <- which(dat$row == 7)
row8 <- which(dat$row == 8)
row9 <- which(dat$row == 9)
row10 <- which(dat$row == 10)
blk_r1 <- c(row1,row2)
blk_r2 <- c(row3,row4)
blk_r3 <- c(row5,row6)
blk_r4 <- c(row7,row8)
blk_r5 <- c(row9,row10)

dat$blk_r[blk_r1] <- 1
dat$blk_r[blk_r2] <- 2
dat$blk_r[blk_r3] <- 3
dat$blk_r[blk_r4] <- 4
dat$blk_r[blk_r5] <- 5

dat$blk_r <- as.factor(dat$blk_r)

#blks
dat <- data.frame(dat, blks = 1)

blk1 <- grep("[1-5]{1}_[1-2]{1}", dat$spindex) #grep returns the row indices that match the regular expression
blk2 <- grep("[1-5]{1}_[3-4]{1}", dat$spindex)
blk3 <- grep("[1-5]{1}_[5-6]{1}", dat$spindex)
blk4 <- grep("[1-5]{1}_[7-8]{1}", dat$spindex)
blk5 <- grep("[1-5]{1}_[9,10]{1}", dat$spindex)
blk6 <- grep("[6-9,10]{1}_[1-2]{1}", dat$spindex)
blk7 <- grep("[6-9,10]{1}_[3-4]{1}", dat$spindex)
blk8 <- grep("[6-9,10]{1}_[5-6]{1}", dat$spindex)
blk9 <- grep("[6-9,10]{1}_[7-8]{1}", dat$spindex)
blk10 <- grep("[6-9,10]{1}_[9,10]{1}", dat$spindex)

dat$blks[blk1] <- 1 #here we apply the recoded value "1" to the matching row indices
dat$blks[blk2] <- 2
dat$blks[blk3] <- 3
dat$blks[blk4] <- 4
dat$blks[blk5] <- 5
dat$blks[blk6] <- 6
dat$blks[blk7] <- 7
dat$blks[blk8] <- 8
dat$blks[blk9] <- 9
dat$blks[blk10] <- 10

dat$blks <- as.factor(dat$blks)

Removing NAs

x <- data.frame(a = 1:10, b=c(NA,1:5,NA,1:3))
na.exclude(x)

Making histograms

hist()
#assigning specific break points
brks = c(0.1,0.2,0.3)
hist(x, breaks = brks)
#however R will only plot if there are points in the cells and will only plot over the range of points.  this
#makes it difficult to make histograms with consistent scales.  In these cases I use barplot() instead.

Making barplot() with Standard Error of the Mean bars (SE bars)

means <- tapply(dat$fin.gall.dia, list(dat$trmt, dat$fin.morph), mean, na.rm = T)
barplot(means[c(1,2),c(1:3)], beside = T)

means <- tapply(dat$fin.gall.dia, list(dat$trmt, dat$fin.morph), mean, na.rm = T)
means
labs <- list(c("Young","Mature"),c("Cushion", "Flat", "Irregular", "Modesta"))
dimnames(means) <- labs

s.size <- tapply(dat$fin.gall.dia, list(dat$trmt, dat$fin.morph), length)
sr.ss <- sqrt(s.size)
se <- means/sr.ss
se

#calculate probabilities of pupae or adults on Apex or Mature
p.ad.d <- tapply(dat$fin.m.stg, list(dat$trmt, dat$fin.morph), length)
p.ad.n <- tapply(which(dat$fin.m.stg == "p" | dat$fin.m.stg == "a"), 
   list(dat$trmt[which(dat$fin.m.stg == "p" | dat$fin.m.stg == "a")], 
   dat$fin.morph[which(dat$fin.m.stg == "p" | dat$fin.m.stg == "a")]), length)
p.ad <- round(p.ad.n/p.ad.d*100, digits = 0)
p.ad #percent surviving to pupae or adults

barplot(means[c(1,2),c(1:3)], horiz = F, beside = T, ylab = "Mean Final Gall Diameter (mm)",
 ylim = c(0,5), legend.text = T, col = c("grey", "white"))
#add standard error of the mean bars (i.e., +/- x/sqrt(n))
arrows(1.5,means[1] - se[1], 1.5, means[1] + se[1], angle=90, code=3, length=0.03)
arrows(2.5,means[2] - se[2], 2.5, means[2] + se[2], angle=90, code=3, length=0.03)
arrows(4.5,means[3] - se[3], 4.5, means[3] + se[3], angle=90, code=3, length=0.03)
arrows(5.5,means[4] - se[4], 5.5, means[4] + se[4], angle=90, code=3, length=0.03)
arrows(7.5,means[5] - se[5], 7.5, means[5] + se[5], angle=90, code=3, length=0.03)
arrows(8.5,means[6] - se[6], 8.5, means[6] + se[6], angle=90, code=3, length=0.03)
text(1.5,means[1] + se[1] + 0.15, p.ad[1])
text(2.5,means[2] + se[2] + 0.15, p.ad[2])
text(4.5,means[3] + se[3] + 0.15, p.ad[3])
text(5.5,means[4] + se[4] + 0.15, p.ad[4])
text(7.5,means[5] + se[5] + 0.15, p.ad[5])
text(8.5,means[6] + se[6] + 0.15, p.ad[6])

Plotting multiple graphs in one device frame

par(mfrow = c(2,2)) #prints four graphs (by row) in a 2 x 2 space
par(mfcol = c(2,2)) #prints four graphs (by column) in a 2 x 2 space

# For finer control over widths and arrangement use layout #
# One figure in row 1 and two figures in row 2
# row 1 is 1/3 the height of row 2
# column 2 is 1/4 the width of the column 1
attach(mtcars)
layout(matrix(c(1,1,2,3), nrow = 2, ncol = 2, byrow = TRUE),
   widths=c(3,1), heights=c(1,2))
hist(wt)
hist(mpg)
hist(disp)

#to practice getting it right before plotting
nf <- layout(matrix(c(2,0,1,3),nrow = 2, ncol = 2,byrow=TRUE), c(3,1), c(1,3), TRUE)
layout.show(nf)

Expand a dataset, also see expand() function

#say you have a dataset that has a count for observations, but all the other variables have the same
#value. For example like this:

> dat <- data.frame(count = c(2,2,3), color = c("green","blue","orange"), fecund = c(8,9,10),stringsAsFactors = F)
> dat
  count  color fecund
1     2  green      8
2     2   blue      9
3     3 orange     10

#you want to expand this so that each count is 1 and therefore you can eliminate it.

Method 1:

> rec1 <- dat$count
> rec <- rec1
> rec <- ifelse(rec1 == 0, rec <- 1,rec <- rec1) #need to set zeros to 1 so they are not lost
 
> dat.ex <- data.frame(count = 1:sum(rec), color = "green", fecund = NA, stringsAsFactors = F) #new dataframe to fill

> j <- 0
> for (i in 1:length(dat$count)){
+  rp <- rec[i]#fill rp rows
+  st <- j + 1 #start filling here
+  dat.ex[st:(j+rp),] <- dat[i,]#fill rows in
+  j <- j+rp #update j
+ }
 
> dat.ex
  count  color fecund
1     2  green      8
2     2  green      8
3     2   blue      9
4     2   blue      9
5     3 orange     10
6     3 orange     10
7     3 orange     10

Method 2: 
# convert data to zeros and ones
dat1 <- dat[rep(1:nrow(dat), dat$alive), 1:5]
dat2 <- dat[rep(1:nrow(dat), dat$dead), 1:5]
dat1 <- data.frame(dat1, surv = 1)
dat2 <- data.frame(dat2, surv = 0)

dat3 <- rbind(dat1,dat2)

Sample a vector of subjects or treatments without replacement to make randomized complete blocks

# method 1
#this makes 5 unreplicated blocks with each block having each genotype in random order

x <- c("10_12","10_2","10_5","10_7","10_8","10_9","11_1","11_10","11_11","11_4","11_6","11_7","12_2","12_3","20_1",
"20_10","20_3","20_4","20_7","20_8")

a <- sample(x,20,replace=F)
b <- sample(x,20,replace=F)
c <- sample(x,20,replace=F)
d <- sample(x,20,replace=F)
e <- sample(x,20,replace=F)

y <- c(a,b,c,d,e)
y <- matrix(y, nrow = 10, byrow=FALSE)#byrow=FALSE makes it fill by column
y <- data.frame(y)
y

# method 2

dist <- c("N-No30","N-HP30","N-HP24","N-No24","M-No30","M-HP30","M-HP24","M-No24")
block <- data.frame(c1 = NA,c2 = NA,c3 = NA,c4 = NA)
for (i in seq(1,100,by = 2)){ #only use odd numbers so correctly shifts down when filling dataframe "block"
  samp <- sample(dist,8,replace = F)
  block[i,1:4] <- samp[1:4]
  block[i+1,1:4] <- samp[5:8] 
  }
dim(block)
head(block)

write.table(block, file = "randomized_block.txt",
	append = FALSE, quote = FALSE, sep = "\t",row.names = T,
            col.names = T)

Write a file to a folder (i.e., output data to a file)

write.table(y, file = "C:/Documents and Settings/Jeremy/Desktop/Junk/genotypes/2009/random.txt",
	append = FALSE, quote = FALSE, sep = "\t",row.names = FALSE,
            col.names = FALSE)

Capture output to a text file

zz <- file("gradients_v1.txt", open = "wt")
sink(zz)
sink(zz, type = "message")
summary(gam2)
## back to the console
sink(type = "message")
sink()

Specify distance between graphs in a multiframe plot, (i.e., the margins around plots)

#the "mar" argument specifies this with c(bottom,left,top,right)
#within the par() function
par(mfrow=c(2,5), mar = c(5, 2.5, 4, 0.4) + 0.1)

Specify distance of the axis title from the axis tick labels
#first increase the size of the left margin
par(mar = c(5.1,5.1,4.1,2.1))
#then use mtext() to redraw the axis title after specifying ylab="" in the plot function
mtext("axis title text", side = 2, cex = 1.6, ...)

Specify distance the labels are from the axis
par(mgp = c(3,1,0)) #3 is the axis label, 1 is the tick labels, 0 not sure main title maybe

Suppressing tick labels, no axis tick labels
plot(rnorm(20), axes = FALSE)
 axis(1, labels = FALSE)
 axis(2, labels = FALSE)
 box(bty = "l") #also A character string which determined the type of box which is drawn about plots. If bty is one of "o" (the default), "l", "7", "c", "u", or "]" the resulting box resembles the corresponding upper case letter. A value of "n" suppresses the box.

Suppressing top and right box lines on a graph (or suppress right axis)
bty = "l" #suppresses top and right lines of the bounding box
bty = "n" #suppresses whole bounding box and leaves just the axes x and y
#both these are passed in the plot() function to par()

Putting tick labels on the inside or outside of the axis
tck = 0.01 #positive length value goes on the inside
tck = -0.01 #negative length value goes on the outside
#both these are passed in the plot() function to par()

Attributes are in lists and character format to convert to numeric vectors:

r_names <- as.numeric(unlist(attributes(residuals(fit3))))

Sorting a dataframe in ascending order

crdat <- crdat[order(crdat$g, crdat$blks), ]

Using tapply to get means of groups
means <- tapply(dat$nM, list(dat$gld_prs), mean)

To get means, standard deviations (sd), and coefficient of variation (cv), of two columns by row, as a vector
dat <- data.frame(dat, abs_m_pp = rowMeans(dat[,20:21]))
dat <- data.frame(dat, abs_sd_pp = apply(dat[,20:21],1,sd)) #the 1 is by row, 2 is by column
dat <- data.frame(dat, abs_cv_pp = (dat$abs_sd_pp/dat$abs_m_pp)*100)

Comparing means with Tukey's Honestly Significant Difference Test
aov.1 <- aov(dat$achenes ~ dat$geno)
summary(aov.1)
TukeyHSD(aov.1)
plot(TukeyHSD(aov.1))

Here's a photo of the pch attribute of the plot function (i.e., graphic symbols)
Note for filled symbols:  use bg = "red" to change the fill color of pch = 21:25

Specify all two-way (2-way) or greater interactions in an lm model
lm1 <- lm(cap ~ (pesticide + row + col)^2, data)

Making a legend for your plot symbols (cex for letters, pt.cex for symbol size)
legend(2.6, 0.33, c("flat", "cushion", "irregular", "crescent"),
  cex = 1.2, lty = c(1,2,1,2), lwd = c(1.7, 2, 2, 2),
  pch = c(1, 1, 16, 17), bty = "n", col = c("#9E0508", "black", "grey", "orange"),
  pt.cex = c(2,2,3,1.5))

Making a repeating vector (i.e., repeat a number in a vector):
a <- rep(1, 10) #makes a vector of 1's of length ten.

Ordination techniques in R: Principal components analysis; Non metric multidimensional scaling; Discriminant function analysis

# this is a good site too for PCA:
# http://www.r-bloggers.com/computing-and-visualizing-pca-in-r/

##### get data
r.file <- "L:/JJHeath/TLC_ag/To_Score/data_hand/blue_i.txt"
dat1 <- read.delim(r.file,header=T)

#make a matrix to track the names
datm <- as.matrix(dat1[,35:49])
dimnames(datm)[[1]] <- as.character(dat1$con)
datm <- data.frame(datm)

pcdat <- princomp(datm) #principal components analysis (PCA); use covariance matrix
#pcdat <- princomp(datm,cor=TRUE) #use correlation matrix

summary(pcdat) #print standard deviation and proportion of variances for each component
screeplot(pcdat) #It will plot screeplt 
loadings(pcdat) #information how much each variable contribute
# to each component. For principal components you can ignore loading
# subsection of the output from this command

#make vectors for color and pitch
color <- as.character(rep("g",35))
for (i in 1:35){
 if (dat1$morph[i] == "f") {color[i] <- "#9E0508"}
 else if (dat1$morph[i] == "i") {color[i] <- "grey"}
 else if (dat1$morph[i] == "cr") {color[i] <- "orange"}
 else if (dat1$morph[i] == "c") {color[i] <- "black"}
 else {color[i] <- "green"}
 }

library(rgl)
plot3d(pcdat$scores[,1:3], col=color, size=2, type="s")
biplot(pcdat) #or biplot.princomp(pcdat,scale=1) - It will give you biplot

#non-metric multidimensional scaling
library(MASS)
disnme <- dist(datm, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
disnmm <- dist(datm, method = "manhattan", diag = FALSE, upper = FALSE, p = 2)
car.mds <- isoMDS(disnmm)
plot(car.mds$points, type = "n")
text(car.mds$points, labels = as.character(dat1$con))

# Linear Discriminant Analysis with Jacknifed Prediction
library(MASS)
attach(dat1)

fit <- lda(morph ~ p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + p11 + p12 + p13 + p14 + p15,
  data = dat1, na.action="na.omit", CV=F)
fit # show results
plot(fit, col = color, cex = 2)
detach(dat1)

Ordination with FactoMineR
#try FactoMineR package for PCA analysis
library(FactoMineR)
#below need geno_ID in the data argument, then call it quali.sup
res.pca = PCA(datpc[,1:62], scale.unit=TRUE, ncp=15, quali.sup=1, graph=T) 

concat <- cbind.data.frame(datpc[,1],res.pca$ind$coord)
ellipse.coord <- coord.ellipse(concat, bary = TRUE)

color2 <- c("#9E0508", "grey", "orange", "black", "green", "blue", "red", "purple", "yellow", "brown")

plot.PCA(res.pca, axes=c(1, 2), choix="ind", habillage=1, col.hab = color2, ellipse=ellipse.coord,
 label = list("quali"), cex = 1, cex.axis = 2,
 cex.label = 2,
 cex.lab = 1, ylim = c(-4,2), xlim = c(-3,3))
plotellipses(res.pca, level = 0.99)

##### another using FactoMineR

# I like the FactoMineR package (remember R is case sensitive)
# Here's a good resource:
# https://www.youtube.com/playlist?list=PLnZgp6epRBbTsZEFXi_p6W48HhNyqwxIu&feature=view_all
# so let's load it, you may need to download and install it first
library(FactoMineR)

# using the iris dataset which should already be install with R
# type iris to see it or head(iris) to see the first few rows
dat <- data.frame(Species = iris[,5],iris[,1:4]) #need classification variable in first column
str(dat)
names(dat)

#call the first column quali.sup
res.pca <- PCA(dat[,1:5], scale.unit=TRUE, ncp=15, quali.sup=1, graph=T) 
plotellipses(res.pca, level = 0.99, means = T) # means = T specifies confidence ellipse on the means
plotellipses(res.pca, level = 0.99, means = F) # means = T specifies confidence ellipse on the observations

# if you get the error:
# Error in if (scale[1] > 0) r <- r/scale[1] : 
# missing value where TRUE/FALSE needed

# it is probably because you don't have enough samples in one of the categories
# like n < 2

#### Also, FactoMineR plotellipses() uses the lattice package.  Make sure to use this when making plots or arguments will not work.

Argument settings for the plot() function for the size of text in graphs like axis labels and axis annotations, and line width:
cex.axis = 1.5 #for annotations (i.e., numbers on the axis)
cex.label = 1.5 #for x axis in category plots
cex.lab = 1.5 #this is for the axis label (i.e., the name of the axis)
abline(lm1, lwd = 2, lty = 3)# lwd = line width, lty = line type (dotted, solid, etc.)
#las = 2 makes bar labels vertical/perpendicular to the x-axis
#las = 1 makes bar labels horizontal/parallel to the x-axis
boxplot(dat$prod ~ dat$carid*dat$langle, las=2)

Make a return or new line when printing text on a graph or figure
text(30,1.2,"R-squared = 0.64\nP = 0.006", cex = 1.5) #use "\n"

Rotate text when printing on a graph, vertical text
text(30, 1.2, "R-squared = 0.90", srt = 90)# 90 degree rotation

Complex expression as an axis label with superscripts
plot(...,ylab = expression(paste("Intensity x area (x ", 10^5, ")", sep = "")),...
#below makes superscript and inserts the variable R-squared because expression doesn't work with variables/objects
text(1.5,4,bquote(R^2 == .(R2Val)), cex = 2.25)

Change the order of the levels of a factor so they will plot in the order you want. Reorder factor levels.
age2 <- factor(dat$age, levels=c("Y","M","NT"))

Delete a level of a factor after you removed it from the data set
dat$morph <- factor(dat$morph)

Plot 3 dimensional surface (3D)
library(mgcv)
library(rgl)
scatter3d(cres_d, prob_platy, mh, fit="linear", residuals=TRUE,
   bg="white", axis.scales=TRUE, grid=TRUE, ellipsoid=FALSE, xlab="cres_d", 
  ylab="prob_platy", zlab="mh")

Export graphs in SVG or vector format using different devices
Note: InkScape will import PDFs just fine from R if use workaround (ii) below.
#default parameters:
devSVG(file = "Rplots.svg", width = 10, height = 8, bg = "white", fg = "black", onefile=TRUE, xmlHeader=TRUE)
#usage
require(RSvgDevice)
devSVG()
plot(1:11,(-5:5)^2, type='b', main="Simple Example Plot")
dev.off()

#for PDFs just save as PDF in the graphic device window.  But often plot symbols do not come out correctly.
Two suggestions for workarounds:
(i)  produce PostScript and then convert to PDF using something  like ghostscript (e.g., ps2pdf)
(ii)  use an almost-but-not-quite opaque colour, e.g., rgb(0, 0, 0, .99) for the points.  If the points are not fully opaque, the character is not used.
Note: (ii) is very easy to do and I like it the best!
(iib) also write this way:  rgb(0,0,0,254, maxColorValue = 255)

#to color the boxplot outliers so they import into InkScape properly use outcol = rgb(0,0,0,0.99)
boxplot(resid(lmcap) ~ age2*stage2, notch=T, las = 2, ylim = c(-5,5),
  ylab = "Capitula mass residuals", cex.lab = 1.5, cex.axis = 1.5,
  outcol = rgb(0,0,0,0.99))

Logistic regression example
#This is an example of testing the effect of relative gall height at maturity of crescent
#morphotypes of Asteromyia carbonifera on risk of parasitism (binomial response: 1 = parasitized, 0 = not)
#by Platygaster solidaginis

read.f <- "H:/Asteromyia/Experiments/Genotypes/2009/WSU_plot/analysis_wsu2010conf/WSU_plot_09_copy_v1a.txt"
dat <- read.delim(read.f, header=T)
str(dat)

#for just crescents
datCr <- dat[dat$morph2 == "cres",]
str(datCr)

#model: risk = 1/1+e^-z where, z = B0 + B1X1
platy <- glm(datCr$p_platy ~ datCr$rghm, family = binomial(logit))

#test whether any better with rghm in the model
anova(platy, test = "Chisq")

#gets you estimates of coefficients B0 and B1 and whether significantly diff from zero
summary(platy)

#gets the locations in "datCr" of the values used to create model, since some were missing (i.e., NAs)
cr_names_platy <- as.numeric(unlist(attributes(fitted(platy))))

#plots the fitted values and the model curve as well as raw values
plot(datCr$rghm[cr_names_platy], fitted(platy), ylim = c(0,1), xlim = c(0,1), ylab = "Risk of Platy attack", xlab = "Relative gall height at maturity")
	text(0.5,0.5, "DF = 1, 1370\nDeviance = 38.8\nP << 0.001") #"\n produces a return"
	points(datCr$rghm[cr_names_platy], datCr$p_platy[cr_names_platy], col = "orange")
	curve(1/(1 + exp(-1*(-4.3128 + 3.7972 * x))), from = 0, to = 1, add = T, lwd = 2, col = "red")

Combine two datasets on a set of index variables. Match up the IDs and add data to dataset.

#make fake datasets or use your own
dat1 <- data.frame(block = c(rep(1,4),rep(2,4)), row = c(1,1,2,2), col = c(1,2,1,2), res = rnorm(8, mean = 10, sd = 1))
dat2 <- data.frame(blockp = c(rep(1,4),rep(2,4)), rowp = c(1,1,2,2), colp = c(1,2,1,2), sp = dat1$res)
#change one value in the dataset to be added to test if working okay
dat2[1,1] <- 0
#reorder to test if working okay
dat2 <- dat2[order(dat2$colp) ,]

#add variables to dat1, which will be changed in the for loop below
dat1 <- data.frame(dat1, blockp = 666, rowp = 666, colp = 666, sp_p = 666.6)

for(i in 1:length(dat1$block)){
  w <- which(dat2$blockp == dat1$block[i] & dat2$rowp == dat1$row[i] & dat2$colp == dat1$col[i])
  if (length(w)>0) {dat1[i,5:8] <- dat2[w,1:4]}
  else {dat1[i,5:8] <- NA}
}

# Also see the merge function (more efficient)

merge(dat1,dat2, all.x = T, all.y = T)

Multivariate analysis of variance: MANOVA and canonical discriminant analysis (ordination) for discrete independent variables (i.e., class = factor, not integer or number)

read.f <- "E:/Car_Quantification/Car_Quant/analysis/AG_data_v2.txt"
dat <- read.delim(read.f, header = T)
str(dat)
levels(dat$morph)

man1 <- manova(cbind(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17)~
	morph, data=dat)
#different types of test on significance
summary(man1, test = "Hotelling-Lawley")
summary(man1, test = "Roy")
summary(man1, test = "Pillai")
summary(man1, test = "Wilks")

#test which responses (carotenoids) differ among morphs
summary.aov(man1)
#test which morphs are different (this doesn't work for some reason?)
man2 <- manova(cbind(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17)~
	morph, data=dat, subset = morph %in% c("cr", "c"))
summary(man2)

#plot canonical discriminant analysis (mainly for discrete independent variables)
require(candisc)
can1 <-candisc(man1, term="morph", ndim=1)
can2 <- candisc(man1, term="morph", ndim=2)
plot(can1)
plot(can2)

Reshape a dataframe in preparation for MANOVA (i.e., change from long to wide format or change one column to many columns)
d <- data.frame(sub2 = dat$sub2, name = dat$name, amount = dat$amount, inject = dat$inject)
d2 <- reshape(d,v.names = "amount", idvar = "sub2", timevar = "name", direction = "wide")

Check all correlations (i.e., correlation matrix) between response variables in a multivariate analysis
#check correlation among variables (gives the correlation coefficients)
cor(dat[,16:32])

#also a special package that makes nice graphs
library(car)
scatterplot.matrix(~cres+cush+flat+irr,  reg.line=lm, smooth=FALSE, span=0.5, diagonal = 'histogram', data=bcwdata)

Round or decrease the number decimal places in a number you want to print
text(1,2.5, round(tst$residuals[1],digits=2))

Converting character dates to real date format
as.Date(x)
#character input for the following is "01-27-11"
date_4 <- as.Date(date_3, format = "%m-%d-%y")

#character input for the following is in this format: "01-Mar-00"
date_check <- as.Date(c("10-Aug-10","10-Sep-10","10-Oct-10"), format = "%d-%b-%y")
 
#use %y for character input that is only 2-digit year
#use %Y with the "Y" capital for 4-digit year character input

Accessing attributes from tapply objects using the unlist function
ags <- tapply(dat$galls, list(date, dat$morph2), sum, na.rm=T)
mas <- tapply(dat$sp_mass, list(date,dat$morph2), mean, na.rm = T)
den <- ags/mas
dags <- data.frame(den, dt = as.Date(unlist(dimnames(den)[1])))

Making a vector of colors, pch, lwd (line widths), or sizes for plotting points
This method is more efficient than using the which() function
col <- c(rgb(0,0,0,254, maxColorValue = 255), rgb(255,165,0,254, maxColorValue = 255),
  rgb(158,5,8,254, maxColorValue = 255), rgb(190,190,190,254, maxColorValue = 255))[dat2$morph]
#do the same for the plot symbols
pch <- c(1,17,1,16)[dat2$morph]
cex <- c(1.5,1,2.2,1)[dat2$morph]
lwd <- c(2,2,2,2)[dat2$morph]

Finding the R code of a function that is in a package
plot.PCA()
getAnywhere("plot.PCA")

Calculate slopes or any linear coeficient for a whole dataset
#####add variable gall growth rate by fitting a linear model to the first two dates and taking the slope

x <- data.frame(d1 = as.numeric(dat$dt_fnd),#dt_fnd is in real date format
 d2 = as.numeric(dat$dt_2),#dt_2 is in real date format
 g1 = dat$dia_fnd,
 g2 = dat$dia_2)
 
slp <- function(x){
 s1 <- 1.0
 d <- x[1:2]
 g <- x[3:4]
 if (sum(is.na(d)) >= 1 | sum(is.na(g)) >= 1){s1 <- NA}
 else if (sum(is.na(d)) == 0 & sum(is.na(g)) == 0) {
 s <- lm(g ~ d)
 s1 <- s$coef[2]
 }
 s1
}

slope <- apply(x, 1, slp)
dat <- data.frame(dat, ggr = slope)

Path analysis and Structural Equation Modeling in R
##Use Onyx (freeware) in combination with the OpenMx R package.  Onyx creates OpenMx code from
##path digrams that with some modification can be used in OpenMx.

require("OpenMx");
manifests<-c("growth","rhiz","cap","stem","treat","loc")
latents<-c("x2","x3","x4","x5")
model <- mxModel("prelim", 
mxData(observed = cov(dat2), type="cov", numObs = length(dat2$growth)),
type="RAM",
manifestVars = manifests,
latentVars = latents,
mxPath(from="growth",to=c("stem","rhiz","cap"), free=c(TRUE,TRUE,TRUE), value=c(1.9,0.6,0.5) , 
  arrows=1, label=c("growth_TO_stem","growth_TO_rhiz","growth_TO_cap") ),
mxPath(from="stem",to=c("cap","rhiz"), free=c(TRUE,TRUE), value=c(0.3,0.05) , 
  arrows=1, label=c("stem_TO_cap","stem_TO_rhiz") ),
mxPath(from="treat",to=c("growth","stem","rhiz","cap"), free=c(TRUE,TRUE,TRUE,TRUE), value=c(-0.3,-6,-5,0.1) , 
  arrows=1, label=c("treat_TO_growth","treat_TO_stem","treat_TO_rhiz","treat_TO_cap") ),
mxPath(from="loc",to=c("growth","stem","rhiz","cap"),free=TRUE,value=c(1,1,1,1),arrows=1,
  label=c("loc_TO_growth","loc_TO_stem","loc_TO_rhiz","loc_TO_cap") ),
mxPath(from="x2",to=c("growth"), free=c(FALSE), value=c(1.0) , arrows=1, label=c("x2_TO_growth") ),
mxPath(from="x3",to=c("stem"), free=c(FALSE), value=c(1.0) , arrows=1, label=c("x3_TO_stem") ),
mxPath(from="x4",to=c("rhiz"), free=c(FALSE), value=c(1.0) , arrows=1, label=c("x4_TO_rhiz") ),
mxPath(from="x5",to=c("cap"), free=c(FALSE), value=c(1.0) , arrows=1, label=c("x5_TO_cap") ),
mxPath(from="cap",to=c("rhiz"), free=c(TRUE), value=c(0.1) , arrows=2, label=c("COV_cap_rhiz") ),
mxPath(from="loc",to="treat",free=TRUE,value = -0.0318981,arrows=2,label="COV_loc_treat"),
mxPath(from="treat",to="treat", free=FALSE, value=var(dat2$treat) , arrows=2, label=c("VAR_treat") ),
mxPath(from="loc",to="loc",free=FALSE,value=var(dat2$loc),arrows=2,label="VAR_loc"),
mxPath(from="x2",to=c("x2"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x2") ),
mxPath(from="x3",to=c("x3"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x3") ),
mxPath(from="x4",to=c("x4"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x4") ),
mxPath(from="x5",to=c("x5"), free=c(TRUE), value=c(1.0) , arrows=2, label=c("VAR_x5") )
)

#options()$mxOptions # list all mxOptions (global settings)
#model <- mxOption(model, "Major iterations",  value = 5000, reset = F)  # may optimize for speed 
#model@options # see the list of options you set

results1 <- mxRun(model)
summary(results1)
#get the correlation residuals
#fix the correlation between treat and loc to zero
residuals <- cov2cor(covar) - cov2cor(results1$objective@info$expCov)
round(residuals, digits=4)

#now I calculate standardized coefficients
#first load the function written by Ryne Estabrook (20 Oct 2010)
#his function name is standardizeRAM, which is in the file: StandardizeRAM.R
stand1 <- standardizeRAM(results1)
#these all come out the same as multiple regression, he also gives standard errors for these but these are not
#typically reported in path analysis and I think there's some controversy about how to calculate them anyway.

#now I need to figure out p-values
#first get the vectors of coefficients and SEs
pvals <- function(results){
  res <- data.frame(summary(results)$parameters)[,c(1:6)]
  names <- res[which(res$matrix == "A"),1]
  betas <- res[which(res$matrix == "A"),5]
  SEs <- res[which(res$matrix == "A"),6]
  ts <- betas/SEs
  n <- summary(results)$numObs
  p <- function(x,n) {2*(pt(-abs(x),df = n-2))}
  ps <- p(ts,n)
  out <- data.frame(names = names, estimate = betas, SEs = SEs, t = ts, P_value = round(ps, digits = 5))
  print(out)
}
pvalues1 <- pvals(results1)

final_est <- data.frame(names = summary(results1)$parameters[,1],
  EstimateNB = summary(results1)$parameters[,5],
  SEsNB = summary(results1)$parameters[,6],
  
  EstimateB = summary(results2)$parameters[,5],
  SEsB = summary(results2)$parameters[,6],
  
  EstimateF = summary(results3)$parameters[,5],
  SEsF = summary(results3)$parameters[,6])
  
final_stand <- data.frame(names = summary(results1)$parameters[,1],
  StandEstimateNB = stand1[,5],
  StandEstimateB = stand2[,5],
  StandEstimateF = stand3[,5])
  
final_pvalues <- data.frame(names = pvalues1[,1],
  PvalueNB = pvalues1[,5],
  PvalueB = pvalues2[,5],
  PvalueF = pvalues3[,5])

Optimize a log transformation and testing for normality
trans <- rep(NA,100)
for(i in seq(1,100, by = 1)){
trans[i] <- shapiro.test(log(rustAUC + i, base = 10))$statistic
}

plot(trans ~ seq(1,100, by = 1))

shapiro.test(log(rustAUC + 40, base = 10))
hist(log(rustAUC + 40, base = 10))

Plotting lots of linear relationships in a lattice plot
library(lattice)
xyplot(is120 ~ resid(nlv) | geno_ID, data = dat,
  type = c("p","r"))

Boxplots for a lot of different levels of discrete factor

bwplot(resid(Smod3) ~ stage2 | geno_ID, data = dat, notch = T)

#or

#all y-axes scales optimized
par(mfrow = c(5,4), mar = c(2, 2.5, 2, 0.4) + 0.1)
for (i in 1:length(unique(dat$name))) {
  boxplot(dat$amount[dat$name == dat$name[i]] ~ 
    dat$inject[dat$name == dat$name[i]], notch = T, main = dat$name[i], las = 2)
  }

How to cite R packages
citation("mgcv") #package name in quotes to find out how to cite it.

Convert times from excel to R time format
dates <- c("03/25/14")
tfin <- dat$time_f #time stop feeding and start searching
tfin <- paste(dates, tfin)
tfin <- strptime(tfin, "%m/%d/%y %H:%M:%S")

Convert times that are positioned around midnight
# make time data when original format is "21:23 or 0:23 or 1:23)

tin <- as.character(dat$time_in) #time put in field
# indices for those before, at, or after midnight
ind.a <- grep("^[0-5]{1}:", tin, value = F) 
ind.b <- grep("^[0-3]{2}:", tin, value = F)
tin[ind.b] <- paste("07/20/14", tin[ind.b], sep = " ")
tin[ind.a] <- paste("07/21/14", tin[ind.a], sep = " ")
tin <- strptime(tin, "%m/%d/%y %H:%M")
tin.num <- as.numeric(tin)
tin.s <- scale(tin.num)

Plot a small plot as an inset in a larger plot

plot(type1 ~ xis, ylim = c(0,0.1), xlab = "Skewness", ylab = "Probability of type 1 error")
 abline(h = 0.05, col = "red")
#add inset of representative histograms for xi values
  par(fig = c(0.2,0.3,0.67,0.77), new = T, las = 1, cex.axis = 0.7, mar = c(0,0,0,0))
  hist(rsnorm(1000, mean = 10, sd = 1, xi = 0.01), xaxt="n",yaxt="n",xlab="", ylab = "", main = "")
  axis(2, mgp=c(2, .4, 0), tcl = -0.3)
  axis(1, mgp=c(2, 0, 0), tcl = -0.3)
  # add titles for the axes
  mtext("Frequency", side=2, line=1.5, cex=0.7, las=0)
#add inset of representative histograms for xi values
  par(fig = c(0.45,0.55,0.35,0.45), new = T, las = 1, cex.axis = 0.7, mar = c(0,0,0,0))
  hist(rsnorm(1000, mean = 10, sd = 1, xi = 1), xaxt="n",yaxt="n",xlab="", ylab = "", main = "")
  axis(2, mgp=c(2, .4, 0), tcl = -0.3)
  axis(1, mgp=c(2, 0, 0), tcl = -0.3)
  # add titles for the axes
  mtext("Frequency", side=2, line=1.5, cex=0.7, las=0)
#add inset of representative histograms for xi values
  par(fig = c(0.8,0.9,0.55,0.65), new = T, las = 1, cex.axis = 0.7, mar = c(0,0,0,0))
  hist(rsnorm(1000, mean = 10, sd = 1, xi = 2), xaxt="n",yaxt="n",xlab="", ylab = "", main = "")
  axis(2, mgp=c(2, .4, 0), tcl = -0.3)
  axis(1, mgp=c(2, 0, 0), tcl = -0.3)
  # add titles for the axes
  mtext("Frequency", side=2, line=1.5, cex=0.7, las=0)

Save and load the workspace to/from the current directory
getwd()
setwd("E:/myspace/myimage")
save.image()
load(.RData) #load a previously saved workspace

Calculate z-scores within a dataframe by grouping factor
dat2 <- data.frame(dat, zfactor = zfactor)
dat2 <- data.frame(dat2[dat2$name == "MeSalicylate",])
#calculate z-scores for MeSA amounts
dat2$sMeSA <- ave(dat2$amount, dat2$zfactor, FUN = scale) 

Plot three dimensional data surfaces 3D, 3d, wireframe
#make grid to predict from model
marZ11 <- seq(from = min(z11.16Ald, na.rm = T), to = max(z11.16Ald, na.rm = T), length.out = 30)
marZ9 <- seq(from = min(z9.14Ald, na.rm = T), to = max(z9.14Ald, na.rm = T), length.out = 30)
grid <- expand.grid(marZ11, marZ9)
names(grid) <- c("z11.16Ald","z9.14Ald")
pred1 <- predict(gam1, grid, type = "response")

# NOTE: when using predict() make sure that the model you use (gam1 in this case) was run with the data = ... argument specified.  If not you'll get an error "Error: variables ... were specified with different types from the fit"

#plot predictions
library(lattice)
wireframe(pred1 ~ grid$z11.16Ald*grid$z9.14Ald)
wf <- wireframe(pred1 ~ grid$z11.16Ald*grid$z9.14Ald, zlim = c(0,1),
  main = "", colorkey = T, drape = T, scales = list(arrows=FALSE),
  #adjust color regions on the surface
  #at=seq(from=0,to=1,length=100), col.regions = col.reg.2,
  ylab = list("Z9-16:Ald", rot = -38), 
  xlab = list("Z11-16:Ald", rot = 34),
  zlab = list("Probability of mating",rot = 94),
  screen = list(z = 40, x = -60, y = 0)#best start z = 40, x = -60, y = 0
  )
print(wf)

Getting the number of combinations/permutations is easy with expand.grid
grd <- expand.grid(one = c("c","t"), two = c("a","g"), three = c("a","g"), four = c("c","t"), five = c("a","t"))
grd$NewCol <- do.call(paste, c(grd[c("one", "two", "three", "four", "five")], sep = ""))
grd
length(unique(grd$NewCol))

Counting up your sample size by a factor without including NAs
tapply(dat$Hv.males, list(dat$trap),
   function(x) length(which(!is.na(x))))

Using the ifelse statement to replace values
#here the values in cors are replaced with NA if they are less than or #equal to 0.5 and left the same otherwise.
ifelse(cors > 0.5, cors, NA)

Making contour plots with the base function contour() and adding marginal histograms or kernel density plots on the margins.

## get needed library

require(mgcv)

## make some foo traits and a fitness function

set.seed(40)
trait.x <- rnorm(mean = 0, sd = 1, n = 1000)
trait.y <- rnorm(mean = 0, sd = 1, n = 1000)
trait.a <- rnorm(mean = 1, sd = 1, n = 100)

## make disruptive selection

u <- 0.3
G <- 0.4

W <- rbinom(1000, 1, exp(u + G*trait.x^2 + G*trait.y^2)/(1 + exp(u + G*trait.x^2 + G*trait.y^2)))
dat <- as.data.frame(list(W = W,tx = trait.x, ty = trait.y))
names(dat)
head(dat)

## run GAM model of fitness function

ff.1 <- gam(W ~ te(tx,ty, k = 10), data = dat, family = binomial)
summary(ff.1)
plot(ff.1)

## Prepare a grid for contour plot

x.mar <- seq(min(trait.x), max(trait.x), length.out = 40)
y.mar <- seq(min(trait.y), max(trait.y), length.out = 40)
grid.1 <- expand.grid(tx = x.mar, ty = y.mar)
grid.1[, "fit"] <- predict(ff.1, newdata = grid.1, type = "response", se.fit = T)
head(grid.1)

## Method 1: Use contour() in conjunction with layout().

# NOTE: Make sure x,y are correct! These are not the x,y values
# generated from expand.grid(), but the sequences fed to expand.grid(),
# which in this case are only length = 40 each. z is a matrix with
# increasing columns corresponding to increasing y-values and rows
# corresponding to increasing x-values. If x, y range is submitted to
# expand.grid in that order then the resulting predictions will be in
# the correct order if you fill the z matrix by column (the default).
# The matrix z will correspond in position with the contour plot if it
# is turned 90 degrees counter clockwise.

lay.1 <- layout(matrix(c(rep(c(2,1,1,1),3),0,3,3,3), nrow = 4, ncol = 4, byrow = T))
layout.show(lay.1)

par(mar = c(3,3,3,0.5))
contour(x = x.mar, y = y.mar, z = matrix(grid.1$fit, nrow = 40, ncol = 40), ylim = c(-3.5,3.5), xlim = c(-3.5,3.5))
  abline(h = 0, v = 0)
  points(0.5,0.4, cex = 3, col = "orange", pch = 19)
den.y <- density(trait.y)
den.x <- density(trait.x)
den.a <- density(trait.a)
#reverse x and y and xlim to rotate plot 
plot(den.y$y, den.y$x, type='l',ylim=c(-3.5,3.5), xlim = c(1,0)) 
  points(den.a$y, den.a$x, type='l',ylim=c(-3.5,3.5), xlim = c(1,0), col = "red")
#change margins so plots line up and are same size
par(mar = c(3,3,0.5,0.5))
#reverse ylim to plot the line down instead of up
plot(den.x$x, den.x$y, type='l',xlim=c(-3.5,3.5), ylim = c(1,0)) 

Assign a vector of character names to a list of objects. Allows to convert vector of strings to objects.

assign("para.ns", mget(paste(para.names, ".n", sep = "")))

Run a script within a script

setwd("...")
source("...")#R script you want to run

Plot on two different y axes. First on one axis then the other.

par(mar = c(5,5,2,5))
with(d, plot(x, logp, type="l", col="red3", 
             ylab=expression(-log[10](italic(p))),
             ylim=c(0,3)))

par(new = T)
with(d, plot(x, n, pch=16, axes=F, xlab=NA, ylab=NA, cex=1.2))
axis(side = 4)
mtext(side = 4, line = 3, 'Number genes selected')
legend("topleft",
       legend=c(expression(-log[10](italic(p))), "N genes"),
       lty=c(1,0), pch=c(NA, 16), col=c("red3", "black"))