Homework 2:

For this assignment you will practice

Calculating summary statistics
Extrating information from a list
Making plots

In this document I

Load the rat bladder data
Make table 30.1
Make Figure 30.1
Make the top of Table 30.2

Your assignment will be

Calculate the other information in Table 30.2 (except “F TEST TO COMPARE VARIANCES”)
Make Figure 30.3
Put all the necessary information into a .Rmd file so that you can “knit” together a Word file.

The Data

Data source

Original Data source: Frazier et al 2006. Effects of gender, age and hypertension on beta-adrenergic receptor function in rat urinary bladder. Naunyn-Schmiedeberg’s Archives of Pharmacology. 373: 300-309

Used in: Motulsky 2nd Ed, Chapter 30, page 220, Table 30.1. Maximal relaxaction of muscle strips of old and young rat bladders stimualted w/ high concentrations of nonrepinephrine (Frazier et al 2006). Response variable is %E.max

Load raw data

“Stack” raw data

#make stacked dataframe
library(reshape2)
df.rat.stack <- melt(df.rats,
                     variable.name = "rat.age",
                     value.name = "percent.E.max")

## No id variables; using all as measure variables

Data tables

Basic data table

A basic table using pander::pander()

Table 30.1. Maximal relaxation of muscle strips of old and young rat bladders with high concentrations of norepinephrine (Frasier et al., 2006).

## Warning: package 'pander' was built under R version 3.3.2

old	young
20.8	45.5
2.8	55
50	60.7
33.3	61.5
29.4	61.1
38.9	65.5
29.4	42.9
52.6	37.5
14.3	NA

Figures

Note: use “stacked” data made with melt() function.

ylim = adjust the y axis
col = sets colors; since there are 2 groups there are two values within the c(…)
pch = sets the point type; 16 happens to be a filled circle
labels = in beeswarm() sets the group lables
ylab and xlab set the axis labels
mext() is used to put a label closer to the y axis
summarBy() is used to calcualte means
points() is used to overlay the means on the raw data
pch = 3 sets the plotting symbol to be a “+”
cex = is “character expansion” and increase the size of the point
we could reproduce the plot more closely with a vertical line for the mean, but this would take more code

library(beeswarm)

## Warning: package 'beeswarm' was built under R version 3.3.2

par(mar = c(3,3.1,1,1))
beeswarm(percent.E.max ~ rat.age, #formula
         data = df.rat.stack,
         ylim = c(0,75), #set y axis limits
         col = c(1,4), #colors; 1 = black, 4 = blue
         pch = 16,      # plotting symbol
         labels = c("Old","Young"),
         xlab = "",     #make xlab blank
         ylab = "")     #makle ylab blank

mtext(text = "%Emax",
      side = 2,
      line = 2)

# calcualte means
library(doBy)
means <- summaryBy(percent.E.max ~ rat.age, #formula
         data = df.rat.stack,
         na.rm = T)

points(percent.E.max.mean ~ rat.age, 
       data = means, pch = 3, cex = 4)

Plotting figures side by side with par()

par(mfrow = c(1,2)) sets up 2 plots, side by side

#set up 1 x 2 plotting configuration
par(mfrow = c(1,2))

#First figure
beeswarm(percent.E.max ~ rat.age, #formula
         data = df.rat.stack,
         ylim = c(0,75), #set y axis limits
         col = c(1,4), #colors; 1 = black, 4 = blue
         pch = 16,      # plotting symbol
         labels = c("Old","Young"),
         xlab = "",     #make xlab blank
         ylab = "")   

#Second figure
beeswarm(percent.E.max ~ rat.age, #formula
         data = df.rat.stack,
         ylim = c(0,75), #set y axis limits
         col = c(3,6), #colors; 1 = black, 4 = blue
         pch = 1,      # plotting symbol
         labels = c("Old","Young"),
         xlab = "",     #make xlab blank
         ylab = "")

t-test objects

I can reproduce Table 30.2 by extracting information from a t-test object.

Note: we are assuming equal variances using var.equal = T

rat.t.test <- t.test(percent.E.max ~ rat.age, #formula
         data = df.rat.stack,
         var.equal = T)


#Extract normal t-test stuff
rat.p <- rat.t.test$p.value

#round p
rat.p <- round(rat.p,4)

t <- rat.t.test$statistic
df <- rat.t.test$parameter

#round t
t <- round(t,3)

#put some text into objects
rat.stars <- c("**")
tails <- "Two tailed"

#Put df and t into a vector together using paste()
df.t <- paste(t,df, sep = ", ")

Make Table 30.2

Don’t worry about these details

#The names
my.rownames <- c("P value",
              "P value summary",
              "Are means significantly different (P<0.05)",
              "One- or two-tail P value?",
              "t and df")

#the info
statoutput <- c(rat.p,rat.stars,"Yes",tails,df.t)

#put into dataframe
table30.2 <- data.frame(my.rownames,statoutput)


#change names of table
names(table30.2) <- c("UNPAIRED t TEST","")

Plot table

pander(table30.2,justify ="left")

UNPAIRED t TEST
P value	0.003
P value summary	**
Are means significantly different (P<0.05)	Yes
One- or two-tail P value?	Two tailed
t and df	-3.531, 15

What to do for the t-test part of the assignment

All you need to do is calcualte the rest of the things in Table 30.2, ignoring the “F TEST TO COMPARE VARIANCE”.

For the line “Mean +/- SEM of column A” you would jsut do this:

mean.young <- rat.t.test$estimate[2]
n.young <- 8
SD.young <- sd(df.rats$young, 
               na.rm = T) #NOTE: this is key

SE.young <- SD.young/sqrt(n.young)

NOTE 0: There is an NA value among the young rats. This makes R’s mathematical functions mad. You have to tell R to skip it by including “na.rm = T”.

NOTE 1: “SEM” = “standard error of the mean” = standard error. Also note that since the standard error of the raw data doens’t figure directly into the calculations of the t test, its not part of the t-test output and has to be calcualted seperately.

NOTE 2: Motulksy has a typo in Table 30.2 He writes “Column A” but reports the data for Young rats, which are in the 2nd columns of Table of 30.2, which most peole would consider column B.

NOTE 3: To get R^2 you will have to use the lm() function. This requires first saving an lm() object, then saving the summary of that lm object, then accessing the r.squared part of that list. In the example below I access the R squared value for a different dataset as an example.

Example of getting R^2 value

#Load the cats data
library(MASS)
data(cats)

#Run an lm and store the output
cat.lm <- lm(Bwt ~ Sex, data = cats) #formula

#Run a sumamry on the lm output and save it
cat.lm.summ <- summary(cat.lm)
    
#Look at the structure of the summary object
str(cat.lm.summ)

## List of 11
##  $ call         : language lm(formula = Bwt ~ Sex, data = cats)
##  $ terms        :Classes 'terms', 'formula'  language Bwt ~ Sex
##   .. ..- attr(*, "variables")= language list(Bwt, Sex)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "Bwt" "Sex"
##   .. .. .. ..$ : chr "Sex"
##   .. ..- attr(*, "term.labels")= chr "Sex"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(Bwt, Sex)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
##   .. .. ..- attr(*, "names")= chr [1:2] "Bwt" "Sex"
##  $ residuals    : Named num [1:144] -0.36 -0.36 -0.36 -0.26 -0.26 ...
##   ..- attr(*, "names")= chr [1:144] "1" "2" "3" "4" ...
##  $ coefficients : num [1:2, 1:4] 2.3596 0.5404 0.0605 0.0737 38.9975 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "(Intercept)" "SexM"
##   .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
##  $ aliased      : Named logi [1:2] FALSE FALSE
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "SexM"
##  $ sigma        : num 0.415
##  $ df           : int [1:3] 2 142 2
##  $ r.squared    : num 0.275
##  $ adj.r.squared: num 0.269
##  $ fstatistic   : Named num [1:3] 53.7 1 142
##   ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
##  $ cov.unscaled : num [1:2, 1:2] 0.0213 -0.0213 -0.0213 0.0316
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "(Intercept)" "SexM"
##   .. ..$ : chr [1:2] "(Intercept)" "SexM"
##  - attr(*, "class")= chr "summary.lm"

#get the R^2 value, r.squared
cat.lm.summ$r.squared

## [1] 0.274543

#save it to an object
R2 <- cat.lm.summ$r.squared


#round it
R2 <- round(cat.lm.summ$r.squared,3)

Motuksly Chapter 30: Re-Making Tables and Charts

brouwern@gmail.com

January 27, 2017