Welcome to Principles of Biology II (BIOL106L)

This manual is used to guide your coding journey as you learn how to use R programming and it assumes that you have spend some time reading through the laboratory manual in preparation for these activities.

In the POBII lab we will be using the lab computers, but we recommend you to also install both R and Rstudio in your personal computer.

You can install and learn about R and Rstudio following the guidelines on the following website:https://RStudio-education.github.io/hopr/starting.html

Intro to R Programming

Setting up R and RStudio

Let’s make sure you have R set up and ready to go! Here is some code to help you through it. You must have R and Rstudio set up before you can move on. However, these programs will be already available to you in the laboratory computers:

If you have not previously done so, go to this webpage and download R: http://cran.r-project.org/

Then go to this webpage and download Rstudio: http://www.rstudio.com/products/RStuido

Once you have done this, open RStudio

Setting up R Markdown

Now we are going to start by getting R markdown all set up. R markdown is a cool and useful tool that allows you to write code, submit your code, and also upload your results all at one time, as a web browser document! It makes it easier to grade, spot mistakes, and allows you to save your submitted assignments. It can take a bit to learn, but I promise it is worth it!

Here are some helpful resources in aiding you to understand R Markdown:

https://posit.co/wp-content/uploads/2022/10/rmarkdown-1.pdf

To set up the packages required for R Markdown to work do the following:

  1. Open an Rscript and type the following code:
if (!require("rmarkdown")) {
  install.packages("rmarkdown", dependencies = TRUE)
  library(rmarkdown)
}

if (!require("knitr")) {
  install.packages("knitr")
  library(knitr)
}

if (!require("tinytex")) {
  install.packages("tinytex")
  library(tinytex)
  tinytex::install_tinytex()  # Install TinyTeX distribution
}
  1. Now click the Run bottom at the top. This should install all the packages needed for R Markdown to work in your computer.

Our First R Markdown script

Let’s open our first R Markdown file:

  1. In RStudio go to “File” at the top of the window. Then click in “New File” and select “R Markdown”. You could also use he icon that looks like a white page with a little “+” on it on the top left hand corner of the program, and select “R Markdown”

  2. This should have opened an window asking for the title of the file, author names and date. You can call it whatever you prefer or just “R Introduction practice”.

  3. Now this should have opened your code editor or R Markdown script file. This file may likely be “untitled” (see the top of the file). Remember that it is important for you to always save your files in your working directory. Please save this file in your working directory by going to “File” -> “Save as”. You can select your working directory and name it “R introduction”.

  4. Before we move to the next step, please verify that this file was saved in your working directory.

  5. In your newly created R Markdown you can see some default text and guidelines in there. Please delete everything after line 6. But make sure you leave the header (lines 1-5) on there!

Fantastic! You are now ready to start writing code in your R Markdown file!

Practicing basic R commands

If you read the lab manual you are now familiarized with the concepts of packages, libraries, variables, vectors and functions. If not, then please revisit the lab manual.

Activity 1: Create a title in R Markdown & regular text

Remember that in R Markdown you will have a section where you will write regular text and a section where you can write your R code called “R chunks”.

When you write regular text in R Markdown there are a series of coding conventions that you can use. The most important one for our class is the use of the number sign “ # “, which to R Markdown it is telling it to start a title. The amount of number signs you add before the text will indicate the size of the header/title. The more # you add the smaller the tittle will be.

- Write the following title in your R Markdown: “Practicing R commands”

- Under the title write a simple sentence: “Today we will practice different R codes”.

  • Now go to the upper part of RStudio and click “Knit”. This is the way you can print the R Markdown report. If you have issues, please ask your instructor.

Activity 2: Creating R chunks

An “R chunk” is where you can write the R code that will be read and processed by R Studio. This chunks work as the console in R. The “R chunks” are delimited by three grave accents ` and parenthesis as follows {r} :

```{r}

    

```

Please note that inside an “R chunk” the number symbol is used to comment out R code. Let’s see how this works. - Let’s add a variable called mytwo <- 2 inside the R Chunk you just created - Let’s make R print the contents of the variable “mytwo”.

#my first variable
mytwo <- 2 #this assigns the number 2 to the variable called "mytwo"
mytwo    # this prints the values of the variable
## [1] 2
  • If inside the newly created “R chunk” we don’t want to print the contents of the variable “mytwo” we can comment out this code by placing the number symbol before it:
#my first variable
mytwo <- 2 #this assigns the number 2 to the variable called "mytwo"
#mytwo    # this prints the values of the variable

-Can you see the difference?, What happened in the second chunk?. Take 5 minutes to understand what just happened

Activity 3: Understanding Variables and Data frames

R doesn’t just work with basic numbers (also called scalars), but also works very easily with vectors and data frames. Data frames are essentially the same as spreadsheets and vectors represent columns. Vectors in R can be created using an R function, c(). (c stands for concatenate). For instance, say that you have fourth measures for the length of growth rate for an eelgrass. You can make a vector that contains these nine values as follows:

growthrate <- c(4.5,5.5,6.6,7.7)

To print the values inside this vector you can just type the vector’s name

growthrate
## [1] 4.5 5.5 6.6 7.7

If you just want to access the first component of this vector you use brackets:

growthrate[1] #the number inside the bracket indicates the specific vector to be accessed
## [1] 4.5
  • How would you access the last component of this vector (a.k.a the value 7.7)?. Please type the code to do so in a new R chunk and describe what the code is doing

Data frames A data frame is just a combination of multiple vectors. Think of it as having an excel sheet with multiple columns where each column represents a specific vector. Fors instance, let’s create a new vector that can complement the growthrate vector and then combine them into a data frame.

Creating a separate vector called “treename”

treename <- c("azul", "dorado", "verde", "rojo") #creating a separate vector

Combining the vectors treename and growthrate to form a data frame named “mytreedata”:

mytreedata <- data.frame(treename, growthrate)
#checkinng the data frame:
mytreedata
##   treename growthrate
## 1     azul        4.5
## 2   dorado        5.5
## 3    verde        6.6
## 4     rojo        7.7

We can separate pieces of the data frame if we would like. Most commonly, you will want to isolate columns of a data frame to call for specific variables inside your datasets. You can do this with a ‘$’ operator.For instance, let’s access the growth rate of the trees:

mytreedata$growthrate
## [1] 4.5 5.5 6.6 7.7

Now let’s calculate the mean growth rate of the trees:

mean(mytreedata$growthrate)
## [1] 6.075

  • On your own, please add two more numbers to the vector “growthrate” and two more categories to the vector “mytreedata”. Re-create the data frame mytreedata, and re-calculate the mean of growth rate.

  • TAKE a 5 minute break & we will continue with the rest of the activities!

Activity 4: Reading data from excel in R

  1. Download the .csv file from blackboard called tick_scutum.csv. This file has data on the scutum length of ticks sampled from four different cottontail rabbits (1 – 4).

  2. Save the excel file in your working directory. In order to read in your data, the downloaded data MUST be in the same folder in which your R Markdown is saved.

  3. Make sure to set your working directory in RStudio. Go to “Session”, then select “Set Working Directory”, then you can select “Choose directory”.

  4. We can use an R function, read.csv() to read in a spreadsheet that has been saved in a basic text format. The standard is to use comma-separated values files (.csv) where each column is separated by a comma.

  • This function requires two arguments. The first is the name of the file to open. This is also a good time to note that when you are referring to text you will need to surround the text by quotes. The second argument tells the function that there is a header row – a row of column names. We can set this value to T (True). Let’s open the file and save it in a variable called ticks_dat:
ticks_dat <- read.csv("tick_scutum.csv", header = T)
  1. You can view your file with the function “View ()”:
View(ticks_dat)

We often give you data that is in “wide” format. R handles data that is in “long” format more efficiently for the analyses that you do. For more information in wide vs long format visit this website

  1. You can also check whether R is reading your columns correctly. Sometimes R can read a column that is supposed to be “categorical” as “numerical” and vice versa. Always check that R is reading your excel file and columns correctly using the function “str()”:
str(ticks_dat)
## 'data.frame':    37 obs. of  2 variables:
##  $ host         : int  1 1 1 1 1 1 1 1 2 2 ...
##  $ scutum.length: int  380 376 360 368 372 366 374 382 350 356 ...

You can notice that your excel file was read as a data frame with 2 variables or columns: host and scutum.length. Both variables were read as “int” or “integers”. Do you think this is correct?.

  1. Remember that the variable “host” refers the rabbits, and they were called by numbers. However, these numbers reflect a category or ID not a true number to make calculations. Therefore, we need to make sure to tell R to read the variable “host” as categorical or a variable with four factors (or rabbits in this case). We do that with the function as.factor():
ticks_dat$host <- as.factor(ticks_dat$host)
  1. Checking the structure of your data frame:
str(ticks_dat)
## 'data.frame':    37 obs. of  2 variables:
##  $ host         : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 2 2 ...
##  $ scutum.length: int  380 376 360 368 372 366 374 382 350 356 ...

Your data frame now says that it’s variable “host” is composed of four factors. This is more appropriate for performing statistical analyses with this data set in R.

  1. Individual portions of a data frame can also be extracted using matrix notation: dataframename[row,col]. If either the row or column value are missing, it will assume you want all the rows or columns respectively. So, another way to isolate the scutum.length vector (as we did above) would be:
ticks_dat[,2]
##  [1] 380 376 360 368 372 366 374 382 350 356 358 376 338 342 366 350 344 364 354
## [20] 360 362 352 366 372 362 344 342 358 351 348 348 376 344 342 372 374 360
  • Again, if we would like to get the value of scutum lenght for the third row of host:
ticks_dat[3,2]
## [1] 360

  • How would you isolate the scutum length value for the fifth row of host?

  • What would be an easy way to calculate the mean of scutum lenght using this data frame?. Please give it a try on your own

Activity 5: Loading packages & Calculating summary statistics in R

Installing and calling packages

An important part of why we love R so much is because it allows us to visualize our data and perform basic statistical analyses very easy. To that end we need to install some packages and call them into our working space.

The very first time you use a package, get a new computer, or update R, you will need to reinstall your packages to the R program. The general code for installing a package is install.packages(). An example can be seen below. Notice that you do need quotation marks around the package name at this stage. For our statistical analyses we will need a package called Rmisc.

install.packages("Rmisc")

After you install a package, or restart R, you will need to load the library. Note that this must be done each time you open R, not just after installing the package. This is done with the general code library(). Notice you do not need the quotations around the package name at this stage.

library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr

Please note that if you do not install a package and call the package, then the code to calculate the basic statistics will not work. So always remember to install the packages you need and call them with the “library()” statement as described above.

Summary Statistics

There are a number of reasons we produce summary statistics. The first is simply that it gives us a quick view into how our data might look. The second is that when we run statistical tests, they often produce p-values and test statistics that aren’t particularly intuitive. We need descriptive statistics of the data in order to explain the results to someone in a way they will understand.

We should now be ready to estimate basic summary statistics. Let’s say that we are interested in evaluating whether the scutum length of ticks observed varied by rabbit hosts. One of the first things we need to do is check the mean and standard errors of the tick scutum length by rabbit host. If the standard errors of the tick scutum length overlap, this may suggest that the mean tick scutum length did not differ among hosts.

We will use the function summarySE() from the Rmisc package. We need to supply three arguments- the data frame name, the name of the column of data that has the measured values (in quotes) and the name of the column that has your treatment levels (in quotes), usually your categorical variable:

sumstats_scutumlength <- summarySE(data=ticks_dat, groupvars = "host", measurevar = "scutum.length", na.rm=T)
sumstats_scutumlength
##   host  N scutum.length        sd       se        ci
## 1    1  8      372.2500  7.363035 2.603226  6.155651
## 2    2 10      354.4000 11.918240 3.768878  8.525795
## 3    3 13      355.3077  8.919871 2.473927  5.390224
## 4    4  6      361.3333 15.266521 6.232531 16.021231

  • Here we see the various statistics for the four hosts including the mean, standard error (se), standard deviation (sd) and sample size (N). If you would like to have a table like this in your lab report, the easiest way to properly format the table is to copy it into Excel to add cell borders and other required formatting options. Althought, R now have some packages that can do this for you, it is out of the scope of our course.

  • From here, I can compare the specific means given within the sumstats_scutumlength table. In order to calculate the difference between two groups, you subtract value 2 from value 1, and divide by value 2. Try to think about this logically first. If I told you that you got an 85 on the first test, and a 93 on the second test, and I wanted to know how much you improved. You would do (93-85)/85. If you do that in R: (93-85)/85, you get 0.094. This would be equivalent to a 9.4% increase in score, which makes sense.

  • In the case of our tick scutum data, let’s say that we want to evaluate by how much was the mean tick scutum lenght in rabbit host 1 different from rabbit host 2. You would do the following math calculation: (372.2500 - 354.4000)/(354.4000). To get R to do this for you we would need to use the concept we learned above on how to isolate specific values of a data frame. Individual portions of a data frame can also be extracted using matrix notation: dataframename[row,col]:

diffh1_vs_h2 <- ((sumstats_scutumlength[1,3]) - (sumstats_scutumlength[2,3]))/(sumstats_scutumlength[2,3])

diffh1_vs_h2 
## [1] 0.05036682

This value suggest that the tick scutum length in host 1 was 5% higher than in host 2

  • Calculate by how much were the differences in the mean scutum lenght of ticks sampled in host 2 with respect to host 3, and host 3 with respect to host 4

Graded Objective

Okay! We are done for today. Your exit ticket is the pdf of the codes you did today with answers to each of the questions that were asked, along with the wrap-up exercise described below:

  1. Create a title: “Wrap up Exercise”
  2. Download the data from blackboard: “PlantGrowth_example.csv”, save it in your working directory folder and read it in R. This is data of plant weight after exposing the plants to different watering conditions.
  3. Check the structure of the data
  4. Make the variable “water_treatment” as a categorical or factor variable. How many levels are in this variable?
  5. Calculate and print the summary statistics for each water treatment.
  6. Calculate by how much was the control treatment (ctrl) is different from treatment 1 (trt1).

How to submit your work:

You are now ready to save the assignment to a PDF Document! Let me walk you through that. This is how you will be turning in your R projects this semester.

  1. Go to the top of the window and hit the “Knit” button, making sure it is knitting to PDF

  2. Upload your PDF file to blackboard as your submission file.

Intro to R graphing and probability distributions

Data formats

Let’s all remind ourselves of the proper way to format data for analysis and graphing in R!. We often give you data that is in “wide” format. R handles data that is in “long” format more efficiently for the analyses that you do. For more information in wide vs long format visit this website

The figure below shows you the long vs wide format:

Activity 1: Transforming data from wide to long format

  1. The data shown in the table below is reflecting the abundance of fungi collected in different environmental samples across continents. Please identify which are the categorical variables and response variables in this data. Discuss with a partner.
Fungi abundance in environmental samples collected in different continents
Sample America Asia Africa India Europe
Soil 45 12 20 60 34
River water 23 10 25 40 32
Marine water 20 5 40 45 33

  1. Please open excel and type the data shown above in long format into excel. Transform the data into “long” format so that it can be used in R. Think carefully about what your columns should be in this newly formatted file. Talk within your groups to decide what it should be, and we will discuss as a group before you move forward

  2. Save the newly formatted data as a CSV file and upload it to your R environment using the “read.csv” command. Read in the file and name the dataframe “fungi_dat”. If you do not remember how to read a file in R go back to the directions in the previous week.

  3. Remember to check the structure of your data. How is R reading your categorical and response variables?. Please make sure that R is reading your categorical variable as a factor by using the as.factor( ) function.

If your data looks like this:

## 'data.frame':    15 obs. of  3 variables:
##  $ Sample.type    : chr  "Soil" "Soil" "Soil" "Soil" ...
##  $ Continent      : chr  "America" "Asia" "Africa" "India" ...
##  $ Fungi.Abundance: int  45 12 20 60 34 23 10 25 40 32 ...

It should look like shown below after making your variables as factors:

## 'data.frame':    15 obs. of  3 variables:
##  $ Sample.type    : Factor w/ 3 levels "Marine water",..: 3 3 3 3 3 2 2 2 2 2 ...
##  $ Continent      : Factor w/ 5 levels "Africa","America",..: 2 3 1 5 4 2 3 1 5 4 ...
##  $ Fungi.Abundance: int  45 12 20 60 34 23 10 25 40 32 ...

Now your data is ready to use in statistical analysis and for plotting.

Plotting data

Today, we are going to focus on plotting using the package ggplot, which is an incredibly powerful graphics tool. While it can take a bit of work to learn the code, once you learn it, there are very few graphics you can’t make using the same base code.

For more help on the logic of building a graph - click here

This is also helpful - click here

Let’s firts install the package:

#You will need to install the package if you have not already 
install.packages("ggplot2")

Now call the package:

#don't forget to read in the appropriate library!
library(ggplot2)

The first line of code always does the same thing: ggplot(DFname, aes(x=Var, y=Var)). This line uses the command “ggplot” to produce a plot from a specific dataframe (DFname), and you define the X and Y axes of the plot. “aes” means aesthetics - ” deals with the principles of beauty and artistic taste”. When you see “aes”, you are telling ggplot how you want the code to look.

The second line of the code describes the type of plot you are going to make. When you use the code geom_bar() you are telling ggplot that you want it to produce a barplot. When you put geom_bar(stat="summary") you are indicating what the bars will show. The heights of the bars commonly represent one of two things: (1) either a count of cases in each group, or (2) the values in a column of the data frame, such as a representation of mean. By default, geom_bar uses stat=“bin”. This makes the height of each bar equal to the number of cases in each group (#1 above)and will likely produce an error. If you want the heights of the bars to represent values (such as means) in the data, use stat=“summary” and map a value to the y aesthetic. A barplot is one MANY plots you can make with ggplot.

All lines of code following these two are variable depending on what you want the graph to look like. The first two lines are mandatory. You must tell it what data to use, and what type of plot to make. Those lines alone will produce a graph. Everything afterward is to change the appearance of the graph or provide additional information on the plot.

Let’s say that you are interested in evaluating whether there are differences in the abundance of fungi across continents. You would want to plot as your x-axis the “Continents” and as your Y-axis the “Fungi abundance”:

ggplot(fungi_dat, aes(x=Continent, y=Fungi.Abundance)) +
  geom_bar(stat="identity")

This would work perfectly fine if you only wanted to make a barplot with no additional information. However, in this lab, we also want you to include indications of error on your barplots. The best way to add that data to your plot is to use the summarySE function from the Rmisc package we went over a week ago.

Activity 2

  1. Use the summarySE function from the Rmisc package we went over a week ago to calculate the mean fungi abundance found on each continent, and get a measurement of standard error for each continent. Save these results as “stats_continent”. Don’t forget to load the package that runs the summarySE command first!
library(Rmisc)
stats_continent <- summarySE(fungi_dat, measurevar="Fungi.Abundance", groupvars="Continent", na.rm = TRUE)
  1. Now go ahead and view your results! You should have a mean, the standard deviation, the standard error and confidence intervals.
stats_continent
##   Continent N Fungi.Abundance        sd        se        ci
## 1    Africa 3        28.33333 10.408330 6.0092521 25.855725
## 2   America 3        29.33333 13.650397 7.8810603 33.909466
## 3      Asia 3         9.00000  3.605551 2.0816660  8.956686
## 4    Europe 3        33.00000  1.000000 0.5773503  2.484138
## 5     India 3        48.33333 10.408330 6.0092521 25.855725
  1. Now we will use this output to create a new ggplot graphic so that we can incorporate the measurements of error onto our plot. Start by reproducing the above plot using your summarySE results rather than the “fungi_dat” dataframe.
ggplot(stats_continent, aes(x=Continent, y=Fungi.Abundance)) +
  geom_bar(stat="identity")

In your code, the third line is used to add error bars to the plot. We ask you to base your error bars on the standard error. If you remember, you found standard error as part of your summarySE output. Lets take a look at our summary statistics data frame again:

stats_continent
##   Continent N Fungi.Abundance        sd        se        ci
## 1    Africa 3        28.33333 10.408330 6.0092521 25.855725
## 2   America 3        29.33333 13.650397 7.8810603 33.909466
## 3      Asia 3         9.00000  3.605551 2.0816660  8.956686
## 4    Europe 3        33.00000  1.000000 0.5773503  2.484138
## 5     India 3        48.33333 10.408330 6.0092521 25.855725

You can see that I have a column that is labeled “se” and one that is labeled “Fungi.Abundance” and another called “Continent”. Remember to always look at the object you are pulling your data from to ensure that you call the correct variables.

The general code to add errorbars to a graphic is geom_errorbar(aes(ymin=NumVar - se, ymax=NumVar + se), width=0.5). This code tells it to add an errorbar to each bar (each category) that is a certain value below (ymin=NumVar - se) and above (ymax=NumVar + se) the mean for each bar. When you use summarySE, the “se” variable is automatically called that. So you will not have to change it, but you will have to replace “NumVar” with your numerical variable. Think carefully about what variable you want to use here!

If you chose the correct variable, you should have error bars that look like this:

ggplot(stats_continent, aes(x=Continent, y=Fungi.Abundance)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=Fungi.Abundance-se, ymax=Fungi.Abundance+se), width=0.5)

The “width=0.5” part of the command tells R how wide to make the errorbars visually. Try changing it from 0.5 to 0.2 and see what happens!

  1. Up until point you should have a barplot with standard errors. However, you can also modify your X and Y-axes labels to make your figure more professional. Please follow the directions below to change the labels of your axes.

The xlab() and ylab() commands rename your x and y axis to something that makes sense. Because of how we feed data into R, sometimes the way we label variables doesnt make much sense to someone else that isnt familiar with the data. So here, we can rename them! Whatever you put into the parentheses will rename your axis labels, but will not alter the data itself in any way. In this example, I have renamed “Fungi.Abundance” to “Fungi Abundance (%)”. I have also changed the theme of the graphic to change its appearance. You can find more information of ggplot themes here: https://ggplot2.tidyverse.org/reference/ggtheme.html

ggplot(stats_continent, aes(x=Continent, y=Fungi.Abundance)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=Fungi.Abundance-se, ymax=Fungi.Abundance+se), width=0.2) +
  theme_classic() +
  ylab("Fungi Abundance (%)") + 
  xlab("Continents")

If you get stuck with ggplot, GOOGLE IT! There are millions of resources available. It is an incredibly popular tool, and someone before has very likely had the same issue. If you google an error, or a goal that you want to meet, there is likely code online to show you how to do it.

  1. Based on your figure, what can you conclude about your original question? How can you tell whether there were differences in the abundance of fungi among continents?

Graded Objective

  1. Create a new title called “Graded Objective”.

  2. Download the “ChickWeight.csv” file from blackboard. This data is showing the weight of chickens exposed to different diets. Scientists collected the chickens weight in miligrams (mg) at different days (variable=Time).

  3. Check the structure of your data. Think carefully about which variables should be changed to a factor or categorical in R. If so, please make sure do it.

  4. Use this information to get summary statistics of the chicken weight per diet tested.

  5. Create a ggplot barplot using the summary statistics results with error bars at a size of 0.3.

  6. Rename the y-axis to “Weight (mg)”.

  7. Change the theme of your ggplot to anything you’d like.

  8. Please write a paragraph that interprets your graph. Which diet increased the chicken’s weight the most?. Where there differences in the chicken weight among diets?

  9. Print your final plot including all components into your PDF file.

How to submit your work:

  1. You are now ready to save the assignment to a PDF Document! Let me walk you through that. This is how you will be turning in your projects this semester.
  • Go to the top of the window and hit the “Knit” button, making sure it is knitting to pdf. - If you are having trouble knitting as a PDF docuemnt, you can instead knit to HTML. After it produces your HTML webpage, open the file in any browser, then go to File -> Print -> PDF.
  1. Upload your PDF file to blackboard as your submission file.

Hypothesis testing: T-test

In today’s lab we are going to learn how to perform a t-test using R. A t-test (also known as Student’s t-test) is a hypothesis test method that compares the difference in the mean between two measured samples and in the simplest form asks: Are they different?. A t-test may be used to evaluate whether a single group differs from a known value (a one-sample t-test), whether two groups differ from each other (an independent two-sample t-test), or whether there is a significant difference in paired measurements (a paired, or dependent samples t-test). Here we are going to focus on an independent two-sample t-test.

In an independent two-sample t-test your statistical null hypothesis (H0) is that there are NOT significant differences between the means of two samples (H0: μ1 = μ2 ).

We will use as an example data collected from testing an anticancer drug in dogs. In this case scientists discovered a potential anticancer effect from a type of oregano plant only found in the Amazon, so they wanted to test whether the anticancer drug created from this plant is efficient in treating cancer in dogs. They harvested cancerous tissues from different dogs and then exposed them to the drug. They used as a control cancerous tissues that were exposed to just the buffer used to dilute the drug. After a week the scientist then counted the number of cells.

Reading & preparing data

  1. Please download the data dog_anticancerdrug.csv from blackboard.

  2. Make sure to save your data in a working directory.

  3. Open an Rmarkdown file in Rstudio, save the file in your working directory.

  4. Read the “dog_anticancerdrug.csv” into R.

  5. Please check the structure of the data. Is R reading all variables as numerical or categorical?.

## 'data.frame':    60 obs. of  3 variables:
##  $ Dog_ID            : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Treatment         : chr  "drug" "drug" "drug" "drug" ...
##  $ Number_CancerCells: int  9 10 11 12 8 9 13 7 10 9 ...
  1. Based on the description of the data provided above, please answer the following questions:
  1. What is the biological hypothesis that the scientists have?
  2. What is the statistical hypothesis that will help answer the biological hypothesis?
  3. What are the dependent and independent variables in this case? Which variable is categorical?

  1. Make sure the variable that should be categorical is read as a factor. Check your structure again after doing the transformation.

  2. Please calculate the summary statistics.

##   Treatment  N Number_CancerCells       sd        se        ci
## 1   control 30           14.76667 1.546594 0.2823682 0.5775078
## 2      drug 30            9.90000 2.073644 0.3785939 0.7743114

  1. Please estimate by how much was the number of cancer cells different between the control and the drug treatments.

  2. Please write down an interpretation of your result above (in question 9).

Based on the percent difference you calculated above you can only say by how much the control and the drug treatment differed, but you do not know whether this particular percentage is statistically significant. This is when the t-test is needed. Next, let’s first make a plot to visualize our summary statistics and then perform the T-test.

  1. Let’s make a plot using the summary statistics
#calling library ggplot2
library(ggplot2)

myplot <- ggplot(dogdata_sumstats, aes(x=Treatment, y=Number_CancerCells)) +
  geom_bar(stat="identity", color="black") +
  geom_errorbar(aes(ymin=Number_CancerCells-se, ymax=Number_CancerCells+se), width=0.2) +
  theme_classic() +
  ylab("Number of Cancer Cells") +
  xlab("Anticancer drug treatment")

myplot

A prettier graph could be done by modifying some of the components of ggplopt2:

myplot2 <- ggplot(dogdata_sumstats, aes(x=Treatment, y=Number_CancerCells, fill=Treatment)) +
  geom_bar(stat="identity", color="black") +
  scale_fill_manual(values=c("gray", "deepskyblue")) +
  geom_errorbar(aes(ymin=Number_CancerCells-se, ymax=Number_CancerCells+se), width=0.2) +
  theme_classic() +
  ylab("Number of Cancer Cells") +
  xlab("Anticancer drug treatment") +
  theme(axis.title.x = element_text(colour = "black", face = "bold", size=14),
        axis.title.y = element_text(colour = "black", face = "bold", size=14),
        legend.position = "top")
myplot2

You can save your ggplot2 object with the function “ggsave ()” and modify the dimensions:

ggsave(myplot, file="myplot.pdf", width=3.5, height=4)
  1. Based on this plot could you tell whether the differences in the number of cancer cells was statistically different between the two treatments (e.g. control vs drug)?. Please explain your answer.

T-test

Ultimately, you will need to compare the mean number of cancer cells in both groups of dogs to reliably say that the drug that the scientists discovered has anticancer properties. It is in this scenario when you can use a t-test. In this example, our response variable (e.g. number of cancer cells in dogs) is a numerical variable, and the independent variable (e.g. drug treatment) is categorical with two levels: drug administered, drug NOT administered or control group.

In a t-test we are going to: 1) take the difference of our means 2) set an uncertainty we are willing to work with or the threshold to reject a null hypothesis (0.05)
3) show evidence with which we can judge whether the two means are different or not. In our case we are going to return to that 95% number we talked about above. Let’s take the two standard deviations as our rule. If one group’s mean is more than two standard deviations away from another group’s mean we can say they are different (stated mathematically, ∝ = 0.05). So, we interpret our results through our ‘P value’. Our P value is the probability of you getting the number you have (or one more extreme) if the null hypothesis is TRUE. It is similar to the idea of a false positive, which you want to be very small indeed. In this case 5% or less.

Now that your data is loaded, we can perform the t-test to compare the Control and Treatment groups. We will use the function “t.test()”. The code Components for this function are as follows:

  • t_test_result is the name of the data frame storing the output of our test.

  • Number_CancerCells refers to your response variable=the number of cancer of cells counted on a dog, and it is the last column in your data set.

  • Treatment refers to the independent variable= Treatment, with two levels, Control or Treatment. Thie is your Treatment column in your data set.

  • The var.equal = TRUE argument assumes that the variances of the two groups are equal.

  • data = is the name of your data frame, in this case, dogdata

# Perform the two-sample t-test comparing the Control and Treatment groups
t_test_result <- t.test(Number_CancerCells ~ Treatment, data = dogdata, var.equal = TRUE)

Display the results of the t-test

print(t_test_result)
## 
##  Two Sample t-test
## 
## data:  Number_CancerCells by Treatment
## t = 10.304, df = 58, p-value = 1.01e-14
## alternative hypothesis: true difference in means between group control and group drug is not equal to 0
## 95 percent confidence interval:
##  3.921260 5.812073
## sample estimates:
## mean in group control    mean in group drug 
##              14.76667               9.90000

Key Output Values of the t-test:

1. t-value: - The t-statistic (or t-value) represents how far the sample means are from each other relative to the variability in the data.

  • Interpretation: A larger absolute t-value indicates a greater difference between the group means, suggesting a stronger effect and larger effect size. In general, if the absolute t-value is large enough it indicates that the difference between the groups is statistically significant.Example: t = -7 means that the difference between the group means is about 7 times the size of the standard error. This means the difference is likely not due to any error or random chance.

2. Degrees of Freedom (df): This value represents the number of independent values in your dataset that are free to vary. It is used to determine the critical value of the t-distribution.

Interpretation: A higher degrees of freedom generally indicates more data points, which gives a more reliable estimate of the difference between group means.

3. p-value: - The p-value helps you assess whether the difference between the groups is statistically significant.

  • The p-value is determined using the t-distribution and degrees of freedom from a t-table.

Interpretation: - If p-value ≤ 0.05: The difference between the groups is considered statistically significant, meaning there is strong evidence to reject the null hypothesis (which states that there is no difference between the groups).

- If p-value > 0.05: The difference between the groups is not statistically significant, meaning there is insufficient evidence to reject the null hypothesis.

4.Confidence Interval (CI):

  • The confidence interval for the difference in means gives you a range of values within which the true difference between the group means is likely to lie, based on the sample data.

Interpretation: - If the confidence interval does not include 0, it indicates a statistically significant difference between the two groups (because a difference of 0 would mean no difference between the groups).

5. Sample Means: - The mean values of the two groups are often displayed in the output as well. These values represent the average measurement for each group.

Interpretation: The difference in sample means tells you how much higher (or lower) one group is compared to the other.

  1. Please interpret the results of your t-test. Is the difference between control and drug treatments statistically significant?

Graded Objective

  1. Create a new title “Results Section”
  2. Write a paragraph of your results as it is described in the lab manual with the figure and figure caption.

Analysis of Variance (ANOVA)

An ANOVA is used to tests the null hypothesis that the mean of the continuous response variable is the same for every level or factor of the categorical variable (H0: U1 = U2 = U3 …= Um), or that its mean is the same across groups. The alternative hypothesis, is that the means of at least two levels are different from each other. For instance, in the case of our Fastplant Experiment we are interested to test whether different organic fertilizers have the same effect in the plant growth of Brassica rapa plants. In this case the null hypothesis is that the mean growth of B. rapa plants is the same when treated with different types of fertilizers.

Below we will be using a data example, but you should be using your own experimental data.

Class Activity:

  1. Type the data you collected using Table 3 as a template into Excel

  2. Save the Excel file as a “.csv” file:

  • After entering your data, save your spreadsheet as a CSV file:
  • Click on File > Save As.
  • Select CSV (Comma delimited) (*.csv) from the “Save as type” dropdown.
  • Save the file in a folder called “T.Test POB Lab” and name it something like plant_height_data.csv.
  1. Read the file into R and check the structure of your data.
  • Set your working directory in RStudio to the folder where you saved the CSV file. You can do this by navigating to Session > Set Working Directory > Choose Directory.
  • Make sure to save your Rmarkdown file in the same Working Directory folder.You can name it “Anova.rmd”
  • Then, run the following R code to load the data:
#reading the file 
fertdata <- read.csv("fertilizer_data.csv", header=T)
  • Checking structure of the data is important to make sure R is reading your file correctly. You can do this with the function”str()“.
str(fertdata)
## 'data.frame':    90 obs. of  3 variables:
##  $ PlantID  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Treatment: chr  "control" "control" "control" "control" ...
##  $ height   : num  11.4 12.69 12.45 9.82 7.65 ...

Question 1: Can you see which variables are read as categorical and which ones are read as numerical?

  1. Transform your independent variable to be read as a categorical variable in R. You can do this with the function as.factor()
#Making sure that my categorical variable is read as a factor:
fertdata$Treatment <- as.factor(fertdata$Treatment)

#It is always good practice to check your data structure again, to evaluate whether the transformation was effective. 
str(fertdata)
## 'data.frame':    90 obs. of  3 variables:
##  $ PlantID  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Treatment: Factor w/ 3 levels "control","fert1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ height   : num  11.4 12.69 12.45 9.82 7.65 ...

If you have a “Factor” variable with more than three levels something may be wrong with your excel file (check the name of your treatments in the original excel file).

  1. Let’s check how the data looks like to have an idea of what is happening and how is the data distributed by making a boxplot. Remember that a boxplot is a method for demonstrating graphically the locality, spread and skewness of groups of numerical data.
# Load ggplot2 library for visualization
library(ggplot2)

# Create a box plot to visualize the variation in plant height from each fertilizer treatment
ggplot(fertdata, aes(x = Treatment, y = height, fill = Treatment)) +
  geom_boxplot() +
  labs(x = "Fertilizer Treatment", y = "Plant Height (cm)") +
  theme_classic() +
  scale_fill_manual(values = c("Control" = "gray", "fert1" = "deepskyblue", "fert2"="blue"))

  1. Now we are ready to test whether fertilizers have an effect in the plant growth of B. rapa by calculating an ANOVA test. You can do this by using the function aov() :
#calculating ANOVA with the function aov
myanova <- aov(height ~ Treatment, data=fertdata)

To visualize the ANOVA results you need to use the function summary():

#obtaining the ANOVA table with the function summary
summary(myanova)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## Treatment    2  14073    7036    1749 <2e-16 ***
## Residuals   87    350       4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Question 2: What can you conclude from these results?. Please provide a description of your intepretation for your own data

Based on the data we are using in this example the P-value is lower than 0.05, and thus we reject the null hypothesis. Therefore, we would conclude that fertilizer treatment has an effect in plant growth, suggesting that at least one of the fertilizer types administered showed a significantly different effect in plant growth from the other two treatments. However, we do not know from this test which fertilizers had different effects from each other.

Because the ANOVA was statistically significant we will thus perform a post-hoc test. In this case we will use a Tukey HSD Test using the function TukeyHSD(). Tukey honest significant differences test (Tukey HSD) performs multiple t-tests between groups (independent pairwise t-tests), but it controls for the type I error by reducing the cut off value that gives the statistical significance. It considers the ANOVA mean square within calculations to estimate the least amount that means must vary from each other to be significantly different or the Honest Significant Difference (HSD). This test will give you pairwise comparisons between groups. Groupings resulting from the Tukey post-hoc test are usually displayed using letters (“Compact Letter Display”, CLD), where all the groups labeled ‘A’ are not significantly different from each other, ‘B’ are significantly different from ‘A’ but not from other Bs, etc

TukeyHSD(myanova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = height ~ Treatment, data = fertdata)
## 
## $Treatment
##                   diff       lwr      upr p adj
## fert1-control 10.10500  8.870284 11.33972     0
## fert2-control 30.09339 28.858672 31.32811     0
## fert2-fert1   19.98839 18.753670 21.22310     0

Question 3: What can you conclude from the Tukey HSD results?. Make sure to create the groupings and assign letters to each treatment before you proceed to the next step.

  1. Now that we have a clear idea of what treatments may be different from each other, let’s calculate our summary statistics so we can estimate by how much were the treatment means different from each other:
#calling library Rmisc to use the function "summarySE"
library(Rmisc) 

#Calculating summary statistics:
fert_sumstats <- summarySE(data=fertdata, groupvars = "Treatment", measurevar = "height", na.rm=T)
fert_sumstats
##   Treatment  N   height       sd        se        ci
## 1   control 30 10.13452 1.989104 0.3631591 0.7427438
## 2     fert1 30 20.23953 2.102233 0.3838134 0.7849866
## 3     fert2 30 40.22791 1.920937 0.3507134 0.7172895

Obtaining percent differences between means:

control_fert1 <- (fert_sumstats[1,3]- fert_sumstats[2,3])/(fert_sumstats[2,3])
control_fert1
## [1] -0.4992707
control_fert2 <- (fert_sumstats[1,3]- fert_sumstats[3,3])/(fert_sumstats[3,3])
control_fert2
## [1] -0.7480723
fert1_fert2 <- (fert_sumstats[2,3]- fert_sumstats[3,3])/(fert_sumstats[3,3])
fert1_fert2
## [1] -0.4968786

Question 4: Please interpret the percent differences between means. Write three separate statements for each comparison. Which of these comparisons was significant based on the Tukey HSD test you did above?

  1. Laslty, let’s make a plot to have a better visualization of these differences. In the plot we will also add the “Tukey Groupings”
#Let's plot the data based on the summary statistics and adding the Tukey groupings:
myfertiplot <- ggplot(fert_sumstats, aes(x=Treatment, y=height, fill=Treatment)) +
  geom_bar(stat="identity", color="black") +
  scale_fill_manual(values=c("gray", "deepskyblue", "blue")) +
  geom_errorbar(aes(ymin=height-se, ymax=height+se), width=0.2) +
  theme_classic() +
  ylab("Height (cm)") + #remember to include the units you used to measure the plant height!
  xlab("Fertilizer treatment") +
#Tukey group addition: By adding the line of code below to the plot it places text 10 units over the top error # bar (adjust this based on dimensions of Y-axis inplot), and adds CLD text label
  geom_text(label = c("A", "B", "C"), aes(y =(height+se+5),x = Treatment), size = 3) + 
  theme(axis.title.x = element_text(colour = "black", face = "bold", size=14),
        axis.title.y = element_text(colour = "black", face = "bold", size=14),
        legend.position = "top")

myfertiplot

You can save your ggplot2 object with the function “ggsave ()” and modify the dimensions:

ggsave(myfertiplot, file="myfertplot.pdf", width=4, height=4)

Grading Objective

  1. Please create a new title called “Results”
  2. Write a paragraph that describes the results of your analysis
  3. Display your barplot with the Tukey Groupings
  4. Write a Figure caption for the plot