R-markdown
R-markdown
Welcome to week 10! You made it, congratulations. We covered a lot of topics over these last 9 weeks, so we assume you are familiar if not competent with basic R. Today we’re going to put those skills to work and lean R markdown. What is it you ask? R markdown is a system for making dynamic documents. It is text with embedded R code for performing functions like creating summary tables, figures, etc. for reports, posting on webpages. No script for today we’re going all from the webpage. So, buckle up and let’s get started.
Let’s start by creating an R markdown document. First go to File > New File> and choose R markdown.
The following menu will appear
You have the option of choosing your output in HTML, PDF, or Word. For now, let’s choose Word be sure to enter a title and your name. Hit enter and you should see a document with some prepared code. Look closely at the document and you should see that it does not look anything like an R script. There are also strange uses of characters like --- and ```. We see a little R code for example we have summary and plot:
```{r cars}
summary(cars)
```
```{r pressure, echo=FALSE}
plot(pressure)
```
However, they enclosed inside of are highlighted areas. Hmmm… wonder what that means. Well the only way to find out is to click the “Knit” button at the top. It will ask you for the name and location of word file to save and will open the document. First notice that it creates an .Rmd file (aka R Markdown Document) and the corresponding word file. This means that you or anybody else can recreate the analysis using the .Rmd Compare the document with the R code. First you should see the header and text that correspond to the non-shaded area on the R markdown document. You should also see R code and output and a plot with no R code. Now let’s decipher some of that markdown code. It looks like the three dashes delineate the header with key words:
---
title: "My First Markdown"
author: "Jim"
date: "11/28/2022"
output: word_document
---
Go ahead and create another markdown document. This time specify HTML and output and examine. Basically, the same code except that the output is now “html_document”
---
title: "My secondMarkdown"
author: "Jim"
date: "11/28/2022"
output: html_document
---
Go ahead and knit the html document. It looks quite similar with a few things highlighted. I wonder what the output is for PDF? For grins change the html to pdf in the header section and Knit again. What happened? For most of you, you probably got an error because you don’t have LaTex or MiKTeX (MacTeX for macOS) installed. Don’t worry about installing now, this is for future reference. Going back to the html code it looks like the three apostrophe-like things (```) set off the code. They act almost like brackets for functions or for loops. Now let’s examine the text. We see in the document R Markdown is a header and in the code, it is different from other text in that it is preceded by two hashtags and is blue.
## R Markdown
Try removing the hashtags, what happens? It stops being blue and looks like regular text. That is because in R markdown code a hashtag in the text area means header and more hashtags the smaller the header. The table below has the markdown code syntax to do common text functions like bold, subscript, and strike through.
Let’s modify the existing text. This time we’ll create an R markdown document for doing the lesson from last week. In the first text box below the “setup” R code box lets replace the existing text with the following:
## Step 1
We read in the *ground bird data*.
We want the name of the file to be italicized so we add asterisks. Next copy the code for the cars summary and paste immediate below modified code. Delete the “summary(cars)” portion. Next we need to five this piece of code (called a chunk) a title. Click on the settings button at the end of the chunk and a menu will appear as:
Type in the name “Read in data” for the chunk name. You then have several choices as to what this chunk does:
Show output only
Show code and output
Show nothing (run code)
Show nothing (don’t run code)
So there are several options available as well as options for displaying warning messages etc. Choose “Show code and output”. You should now see this in the chunk:
{r Read in data, echo=TRUE}
Cool so the syntax is r “the chunk name” and echo = T means “Show code and output”. Now add the following R code to your chunk:
```{r Read in data, echo=TRUE}
dater<-read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vT9_zoJjILTsVB2bWKNZaDsz7wev0QdIieqpjLheEizSET1-su9gIdcnDLDVYiA9ChjOW4YRjipEzFK/pub?gid=447706920&single=true&output=csv")
```
Below the chunk, lets add the first question to the text box.
##Question 1
Run pairwise correlations between all independent variables (*everything except Density*)
Change the name of the code chunk to “correlation”. Add the code for calculating a correlation. Unlike above, we just want to output the correlations not the code, so echo= F like below.
```{r correlation, echo=FALSE}
cor(dater[,c("Slope", "Elev.m", "Ave.temp", "Ave.ppt", "Pct.graze", "Pct.grass", "Road.dens")])
```
Next question contains a lot of elements a- d. We’re going to skip 3d but will do the rest individually in a separate chunk for each. First task is to plot residuals vs predicted values. Here we will want to just print out the plots so echo= F in the chunk. The text should read:
## Question 3a
For each model *in #2* plot residuals by predicted values.
We then modify the chunk box and add R code for plotting residuals.
```{r plot residuals, echo=FALSE}
# plot residual vs. predicted values
plot(mod1$resid~mod1$fitted, main="Model 1 residual plot", xlab="Fitted values", ylab="Residuals", las=1, pch=15)
plot(mod1$resid~mod2$fitted, main="Model 2 residual plot", xlab="Fitted values", ylab="Residuals", las=1, pch=15)
```
The next question is to calculate AICc for the two models. Here we will be using package MuMIn to extract AICc values so we’ll add a warning to the text and bold the package name. So, in the text box we will add:
## Question 3b
For each model *in #2* plot residuals by predicted values. Note that it requires package **MuMIn**.
Followed by the code chunk:
```{r calculate AICc, echo=FALSE}
AICc1<-MuMIn::AICc(mod1)
AICc2<-MuMIn::AICc(mod2)
D1<-AICc1-min(c(AICc1,AICc2))
D2<-AICc2-min(c(AICc1,AICc2))
df<-data.frame(Model= c("Model 1", "Model 2"), AICc=c(AICc1,AICc2),delta.AICc=c(D1,D2))
df[order(df$AICc),]
```
The third subtask is to create a table with the parameter estimates and 90% confidence intervals, we add that text to the text box as:
## Question 3c
For each model *in #2* create a data frame with parameter estimates and 90% confidence limits.
Below that the code chunk with code and annoying hashtags turned off.
```{r parameters, echo=FALSE, comments = NA}
## grab estimates and SE and combine
parms<-cbind(coef(mod1),sqrt(diag(vcov(mod1))))
## Give the columns the correct names
colnames(parms)<- c("Estimate","Std.Error")
# coerce the matrix into a data frame
parms <- as.data.frame(parms)
# calculate lower and upper 90% CL using t-value 1.64
parms$Lower = parms$Estimate - 1.64*parms$Std.Error
parms$Upper = parms$Estimate + 1.64*parms$Std.Error
#print it out
print("Model 1 parameter estimates with 90% CL")
parms
## grab estimates and SE and combine
parms<-cbind(coef(mod2),sqrt(diag(vcov(mod2))))
## Give the columns the correct names
colnames(parms)<- c("Estimate","Std.Error")
# coerce the matrix into a data frame
parms <- as.data.frame(parms)
# calculate lower and upper 90% CL using t-value 1.64
parms$Lower = parms$Estimate - 1.64*parms$Std.Error
parms$Upper = parms$Estimate + 1.64*parms$Std.Error
#print it out
print("Model 2 parameter estimates with 90% CL")
parms
```
The 4th and 5th questions are related to running logistic regression so we will put the code in a single box and be sure to print the code out. So the test is:
## Question 4 & 5
Fit 2 logistic regression models with the Presence data that have the same independent variables as in #2 above
and the corresponding chunk is lets remove hashtags from output:
```{r logistic regression, echo=T, comment=NA}
dater$Presence<-ifelse(dater$Density>0,1,0)
summary(glm(Presence~ Slope + Elev.m + Pct.graze, data=dater, family='binomial'))
summary(glm(Presence~ Road.dens + Pct.grass + Pct.graze, data=dater, family='binomial'))
```
Now it is time to Knit the markdown and see what we have. The final markdown document can be found here.
This Weeks Assignment
Download the following 2 comma separated (csv) text files for this week's assignments:
coral_count.csv: contains coral counts from 50 quadrats
weather.csv: weather measurements made concurrently with coral counts
Due 1 week from today by 5pm Pacific. Use the 2 data frames (above): coral_count and weather to complete the exact same tasks in week 7 in an R markdown document. You do not have to create a function. You only need to create and output the summary file. Be sure to output code and console output for all steps except the summary file. You should only output the summary table for that task (hint:echo=F). Submit the rmarkdown document to Jim.