Jeromy Anglim's Blog: Psychology and Statistics

Doing Bayesian Data Analysis: A Tutorial with R and BUGS

pdf:
    pandoc talk.md --slide-level 2 -t beamer -o talk.tex
    pdflatex main.tex
    pdflatex main.tex
    -xdg-open main.pdf

\documentclass[t]{beamer}  
\usetheme{Berkeley}                 
\setbeamertemplate{navigation symbols}{}
\title{MY_TALK_TITLE}
\subtitle{MY_TALK_SUBTITLE}
\author{MY_NAME}
\institute{MY_INSTITUTION}
\date{DATE_OF_MY_TALK}

% more preamble...

\begin{document}
\begin{frame}
\titlepage
\end{frame}

\begin{frame}
\frametitle{Outline}
\tableofcontents
\end{frame}

\input{"talk.tex"}

\end{document}

# NAME OF A SECTION
## SLIDE TITLE
* Some point to make
    * Another point
    * Another point
* Some point to make
    * Another point
    * Another point

## SLIDE TITLE
\includegraphics[width=4in]{FIGURE_FILE_NAME.PNG}

function! MarkdownLevel(maxlevel)
    if a:maxlevel >= 1 && getline(v:lnum) =~ '^# .*$'
        return ">1"
    endif
    if a:maxlevel >= 2 && getline(v:lnum) =~ '^## .*$'
        return ">2"
    endif
    if a:maxlevel >= 3 && getline(v:lnum) =~ '^### .*$'
        return ">3"
    endif
    if a:maxlevel >= 4 && getline(v:lnum) =~ '^#### .*$'
        return ">4"
    endif
    if a:maxlevel >= 5 && getline(v:lnum) =~ '^##### .*$'
        return ">5"
    endif
    if a:maxlevel >= 6 && getline(v:lnum) =~ '^###### .*$'
        return ">6"
    endif
    return "=" 
endfunction

au BufEnter *.md  setlocal foldexpr=MarkdownLevel(2)  
au BufEnter *.md  setlocal foldmethod=expr     
au BufEnter *.md  setlocal autoindent

% enables straight single quote
\makeatletter
\let \@sverbatim \@verbatim
\def \@verbatim {\@sverbatim \verbatimplus}
{\catcode`'=13 \gdef \verbatimplus{\catcode`'=13 \chardef '=13 }} 
\makeatother

% enables backticks in verbatim
\makeatletter
{\catcode`\`=13
\xdef\@verbatim{\unexpanded\expandafter{\@verbatim}\chardef\noexpand`=18 }
}
\makeatother

library(xtable)
options(stringsAsFactors = FALSE)
medals <- read.csv("data/medals.csv")
medals$Year <- as.numeric(medals$Year)
medals <- medals[!is.na(medals$Year), ]

# http://www.math.mcmaster.ca/~bolker/emdbook/chap3A.pdf
x <- aggregate(medals$Year, list(Year = medals$Year), length)
names(x) <- c("year", "medals")
x$pos <- seq(x$year)
fit <- nls(medals ~ a * pos^b + c, x, start = list(a = 10, b = 1, 
    c = 50))

plot(medals ~ pos, x,  las = 1, 
        ylab = "Total Medals Awarded", 
        xlab = "Ordinal Position of Olympics",
        main="Total medals awarded 
     by ordinal position of Olympics with
     predicted three parameter power function fit displayed.",
        las = 1,
        bty="l")
lines(x$pos, predict(fit))

medalsByYearByGender <- aggregate(medals$Year, list(Year = medals$Year, 
    Event.gender = medals$Event.gender), length)
medalsByYearByGender <- medalsByYearByGender[medalsByYearByGender$Event.gender != 
    "X", ]
propf <- list()
propf$prop <- medalsByYearByGender[medalsByYearByGender$Event.gender == 
    "W", "x"]/(medalsByYearByGender[medalsByYearByGender$Event.gender == "W", 
    "x"] + medalsByYearByGender[medalsByYearByGender$Event.gender == "M", "x"])
propf$year <- medalsByYearByGender[medalsByYearByGender$Event.gender == 
    "W", "Year"]
propf$propF <- format(round(propf$prop, 2))

propf$table <- with(propf, cbind(year, propF))
colnames(propf$table) <- c("Year", "Prop. Female")

plot(x ~ Year, medalsByYearByGender[medalsByYearByGender$Event.gender == 
    "M", ], ylim = c(0, max(x)), pch = "m", col = "blue", las = 1, ylab = "Total Medals Awarded", 
    bty = "l", main = "Total Medals Won by Gender and Year")
points(medalsByYearByGender[medalsByYearByGender$Event.gender == 
    "W", "Year"], medalsByYearByGender[medalsByYearByGender$Event.gender == 
    "W", "x"], col = "red", pch = "f")

print(xtable(propf$table,
             caption="Proportion of Medals that were awarded to Females by Year"), 
      type="html", 
      caption.placement="top",
      html.table.attributes='align="center"')

cmm <- list()
cmm$medals <- sort(table(medals$NOC), dec = TRUE)
cmm$country <- names(cmm$medals)
cmm$prop <- cmm$medals/sum(cmm$medals)
cmm$propF <- paste(round(cmm$prop * 100, 2), "%", sep = "")

cmm$row1 <- c("Rank", "Country", "Total", "%")
cmm$rank <- seq(cmm$medals)
cmm$include <- 1:10

cmm$table <- with(cmm, rbind(cbind(rank[include], country[include], 
    medals[include], propF[include])))
colnames(cmm$table) <- cmm$row1

print(xtable(cmm$table, caption="Rankings of Medals Won by Country"), 
      "html", include.rownames=FALSE, caption.placement='top',
      html.table.attributes='align="center"')

NOC50Plus <- names(table(medals$NOC)[table(medals$NOC) > 50])
medalsSubset <- medals[medals$NOC %in% NOC50Plus, ]
medalsByMedalByNOC <- prop.table(table(medalsSubset$NOC, medalsSubset$Medal), 
                                 margin = 1)
medalsByMedalByNOC <- medalsByMedalByNOC[order(medalsByMedalByNOC[, "Gold"], 
         decreasing = TRUE), c("Gold", "Silver", "Bronze")]
barplot(round(t(medalsByMedalByNOC), 2), horiz = TRUE, las = 1, 
        col=c("gold", "grey71", "chocolate4"), 
        xlab = "Proportion of Medals",
        main="Proportion of medals won that were gold, silver or bronze.")

listOfYears <- unique(medals$Year)
names(listOfYears) <- unique(medals$Year)
totalNocByYear <- sapply(listOfYears, function(X) length(table(medals[medals$Year == 
    X, "NOC"])))

plot(x = names(totalNocByYear), totalNocByYear, ylim = c(0, max(totalNocByYear)), 
    las = 1, xlab = "Year", main = "Total Number of Countries Winning Medals By Year", 
    ylab = "Total Number of Countries", bty = "l")

ausmedals <- list()
ausmedals$data <- medals[medals$NOC == "AUS", ]
ausmedals$data <- ausmedals$data[, c("Year", "City", "Discipline", 
    "Event", "Medal")]
ausmedals$table <- ausmedals$data

print(xtable(ausmedals$table, 
             caption='List of Australian Medals',
             digits=0),
      type='html', 
      caption.placement='top', 
      include.rownames=FALSE,
      html.table.attributes='align="center"') 

icehockey <- medals[medals$Sport == "Ice Hockey" & medals$Event.gender == 
    "M" & medals$Medal == "Gold", ]
icehockeyf <- medals[medals$Sport == "Ice Hockey" & medals$Event.gender == 
    "W" & medals$Medal == "Gold", ]

# names(table(icehockey$NOC)[table(icehockey$NOC) > 1])

icehockeygs <- medals[medals$Sport == "Ice Hockey" & 
    medals$Event.gender == "M" &
    medals$Medal %in% c("Silver", "Gold"),  c("Year", "Medal", "NOC")]
icetab <- list()
icetab$data <- reshape(icehockeygs, idvar="Year", timevar="Medal",
    direction="wide")
names(icetab$data) <- c("Year", "Gold", "Silver")

print(xtable(icetab$data, 
             caption ="Country Winning Gold and Silver Medals by Year in Mens Ice Hockey", 
             digits=0), 
      type="html",     
      include.rownames=FALSE,
      caption.placement="top",
      html.table.attributes='align="center"')

# if necessary uncomment and install packages.  install.packages('AER')
# install.packages('psych') install.packages('Hmisc')
# install.packages('ggplot2') install.packages('relaimpo')
library(AER)  # interesting datasets
library(psych)  # describe and psych.panels
library(Hmisc)  # describe
library(ggplot2)  # plots: ggplot and qplot
library(relaimpo)  # relative importance in regression

# load the California Schools Dataset and give the dataset a shorter name
data(CASchools)
cas <- CASchools

# Convert grade to numeric

# table(cas$grades)
cas$gradesN <- cas$grades == "KK-08"

# Get the set of numeric variables
v <- setdiff(names(cas), c("district", "school", "county", "grades"))

* district: character. District code.
* school: character. School name.
* county: factor indicating county.
* grades: factor indicating grade span of district.
* students: Total enrollment.
* teachers: Number of teachers.
* calworks: Percent qualifying for CalWorks (income assistance).
* lunch: Percent qualifying for reduced-price lunch.
* computer: Number of computers.
* expenditure: Expenditure per student.
* income: District average income (in USD 1,000).
* english: Percent of English learners.
* read: Average reading score.
* math: Average math score.

str(cas)

## 'data.frame':    420 obs. of  15 variables:
##  $ district   : chr  "75119" "61499" "61549" "61457" ...
##  $ school     : chr  "Sunol Glen Unified" "Manzanita Elementary" "Thermalito Union Elementary" "Golden Feather Union Elementary" ...
##  $ county     : Factor w/ 45 levels "Alameda","Butte",..: 1 2 2 2 2 6 29 11 6 25 ...
##  $ grades     : Factor w/ 2 levels "KK-06","KK-08": 2 2 2 2 2 2 2 2 2 1 ...
##  $ students   : num  195 240 1550 243 1335 ...
##  $ teachers   : num  10.9 11.1 82.9 14 71.5 ...
##  $ calworks   : num  0.51 15.42 55.03 36.48 33.11 ...
##  $ lunch      : num  2.04 47.92 76.32 77.05 78.43 ...
##  $ computer   : num  67 101 169 85 171 25 28 66 35 0 ...
##  $ expenditure: num  6385 5099 5502 7102 5236 ...
##  $ income     : num  22.69 9.82 8.98 8.98 9.08 ...
##  $ english    : num  0 4.58 30 0 13.86 ...
##  $ read       : num  692 660 636 652 642 ...
##  $ math       : num  690 662 651 644 640 ...
##  $ gradesN    : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...

# Hmisc::describe(cas) # For more extensive summary statistics

qplot(expenditure, data = cas) + xlim(0, 8000) + xlab("Money spent per student ($)") + 
    ylab("Count of schools")

round(t(psych::describe(cas$expenditure)), 1)

##            [,1]
## var         1.0
## n         420.0
## mean     5312.4
## sd        633.9
## median   5214.5
## trimmed  5252.9
## mad       487.2
## min      3926.1
## max      7711.5
## range    3785.4
## skew        1.1
## kurtosis    1.9
## se         30.9

# Compute and format set of correlations
corExp <- cor(cas["expenditure"], cas[setdiff(v, "expenditure")])
corExp <- round(t(corExp), 2)
corExp[order(corExp[, 1], decreasing = TRUE), , drop = FALSE]

##          expenditure
## income          0.31
## read            0.22
## math            0.15
## calworks        0.07
## lunch          -0.06
## computer       -0.07
## english        -0.07
## teachers       -0.10
## students       -0.11
## gradesN        -0.17

ggplot(cas, aes(read, math)) + geom_point() + geom_smooth()

# command has strange syntax requiring column numbers rather than variable
# names
partial.r(cas[v], c(which(names(cas[v]) == "read"), which(names(cas[v]) == 
    "math")), which(!names(cas[v]) %in% c("read", "math")))

## partial correlations 
##      read math
## read 1.00 0.72
## math 0.72 1.00

cas$compstud <- cas$computer/cas$students
describe(cas$compstud)

## cas$compstud 
##       n missing  unique    Mean     .05     .10     .25     .50     .75 
##     420       0     412  0.1359 0.05471 0.06654 0.09377 0.12546 0.16447 
##     .90     .95 
## 0.22494 0.24906 
## 
## lowest : 0.00000 0.01455 0.02266 0.02548 0.04167
## highest: 0.32770 0.34359 0.34979 0.35897 0.42083 

qplot(compstud, data = cas)

## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

# Examine correlations between variables
psych::pairs.panels(cas[v])

# (a) students is a count and could be log transformed
cas$studentsLog <- log(cas$students)

# (b) teachers is not the variable of interest:
#   it is the number of students per teacher
cas$studteach <- cas$students /cas$teachers
# (c) computers is not the variable of interest:
#  it is the ratio of computers to students
# table(cas$computer==0) 
# Note some schools have no computers so ratio would be problematic.
# Take percentage of a computer instead
cas$compstud <- cas$computer / cas$students 

# (d) math and reading are correlated highly, reduce to one variable
cas$performance <- as.numeric(
        scale(scale(cas$read) + scale(cas$math)))

m1cor <- cor(cas$performance, cas[c("studentsLog", "studteach", "calworks", 
    "lunch", "compstud", "income", "expenditure", "gradesN")])
t(round(m1cor, 2))

##              [,1]
## studentsLog -0.12
## studteach   -0.23
## calworks    -0.63
## lunch       -0.87
## compstud     0.27
## income       0.71
## expenditure  0.19
## gradesN     -0.16

m1 <- lm(performance ~ studentsLog + studteach + calworks + lunch + 
    compstud + income + expenditure + grades, data = cas)
summary(m1)

## 
## Call:
## lm(formula = performance ~ studentsLog + studteach + calworks + 
##     lunch + compstud + income + expenditure + grades, data = cas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8107 -0.2963 -0.0118  0.2712  1.5662 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.99e-01   4.98e-01    1.80    0.072 .  
## studentsLog -3.83e-02   1.91e-02   -2.01    0.045 *  
## studteach   -1.11e-02   1.59e-02   -0.70    0.487    
## calworks     1.96e-03   2.96e-03    0.66    0.508    
## lunch       -2.65e-02   1.48e-03  -17.97  < 2e-16 ***
## compstud     7.88e-01   3.86e-01    2.04    0.042 *  
## income       2.82e-02   4.89e-03    5.77  1.6e-08 ***
## expenditure  5.87e-05   4.90e-05    1.20    0.232    
## gradesKK-08 -1.21e-01   6.49e-02   -1.87    0.062 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 0.457 on 411 degrees of freedom
## Multiple R-squared: 0.795,   Adjusted R-squared: 0.791 
## F-statistic:  199 on 8 and 411 DF,  p-value: <2e-16 
## 

# calc.relimp from relaimpo package.
(m1relaimpo <- calc.relimp(m1, type = "lmg", rela = TRUE))

## Response variable: performance 
## Total response variance: 1 
## Analysis based on 420 observations 
## 
## 8 Regressors: 
## studentsLog studteach calworks lunch compstud income expenditure grades 
## Proportion of variance explained by model: 79.48%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                  lmg
## studentsLog 0.009973
## studteach   0.016695
## calworks    0.177666
## lunch       0.492866
## compstud    0.025815
## income      0.251769
## expenditure 0.014785
## grades      0.010432
## 
## Average coefficients for different model sizes: 
## 
##                   1X        2Xs        3Xs        4Xs        5Xs
## studentsLog -0.08771 -0.0650133 -0.0558756 -0.0519312 -4.926e-02
## studteach   -0.11918 -0.0861199 -0.0629499 -0.0462155 -3.372e-02
## calworks    -0.05473 -0.0427576 -0.0324658 -0.0233760 -1.535e-02
## lunch       -0.03199 -0.0310310 -0.0301497 -0.0293300 -2.856e-02
## compstud     4.15870  3.0673338  2.2639604  1.6844348  1.287e+00
## income       0.09860  0.0850555  0.0726892  0.0614726  5.140e-02
## expenditure  0.00030  0.0001986  0.0001374  0.0001013  8.061e-05
## grades      -0.45677 -0.3345683 -0.2529014 -0.1981200 -1.628e-01
##                    6Xs        7Xs        8Xs
## studentsLog -4.626e-02 -4.252e-02 -3.833e-02
## studteach   -2.418e-02 -1.687e-02 -1.109e-02
## calworks    -8.399e-03 -2.612e-03  1.962e-03
## lunch       -2.785e-02 -2.718e-02 -2.654e-02
## compstud     1.034e+00  8.828e-01  7.884e-01
## income       4.250e-02  3.477e-02  2.821e-02
## expenditure  6.882e-05  6.206e-05  5.871e-05
## grades      -1.414e-01 -1.291e-01 -1.215e-01

# create a vector of the words that occur in school names
lw <- unlist(strsplit(cas$school, split = " "))

# create a table of the frequency of school names
tlw <- table(lw)

# extract cells of table with count greater than 3
tlw2 <- tlw[tlw > 3]

# sorted in decreasing order
tlw2 <- sort(tlw2, decreasing = TRUE)

# values as proporitions
tlw2p <- round(tlw2/nrow(cas), 3)

# show this in a bar graph
tlw2pdf <- data.frame(word = names(tlw2p), prop = as.numeric(tlw2p), 
    stringsAsFactors = FALSE)
ggplot(tlw2pdf, aes(word, prop)) + geom_bar() + coord_flip()

# make it log counts
ggplot(tlw2pdf, aes(word, log(prop * nrow(cas)))) + geom_bar() + 
    coord_flip()

cas$namelen <- nchar(cas$school)
table(cas$namelen)

## 
## 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 37 38 39 
##  1  4  9 26 28 31 33 27 30 45 38 28 36 30 18 10  5  4  6  3  1  2  2  2  1 

round(cor(cas$namelen, cas[, c("read", "math")]), 2)

##      read math
## [1,] 0.03    0

cas$nameWordCount <- sapply(strsplit(cas$school, " "), length)
table(cas$nameWordCount)

## 
##   2   3   4   5 
## 140 202  72   6 

round(cor(cas$nameWordCount, cas[, c("read", "math")]), 2)

##      read math
## [1,] 0.05 0.01

tlw2p  #recall the list of popular names

## lw
## Elementary      Union       City     Valley      Joint       View 
##      0.983      0.433      0.060      0.040      0.031      0.019 
##   Pleasant        San      Creek        Oak      Santa       Lake 
##      0.017      0.017      0.014      0.014      0.014      0.012 
##   Mountain       Park        Rio      Vista      Grove   Lakeside 
##      0.012      0.012      0.012      0.012      0.010      0.010 
##      South    Unified       West 
##      0.010      0.010      0.010 

# Create a quick and dirty list of popular nature names
naturenames <- c("Valley", "View", "Creek", "Lake", "Mountain", "Park", 
    "Rio", "Vista", "Grove", "Lakeside")

# work out whether the word is in the school name
schsplit <- strsplit(cas$school, " ")
cas$hasNature <- sapply(schsplit, function(X) length(intersect(X, 
    naturenames)) > 0)
round(cor(cas$hasNature, cas[, c("read", "math")]), 2)

##      read math
## [1,] 0.09 0.08

ggplot(cas, aes(hasNature, read)) + geom_boxplot() + geom_jitter(position = position_jitter(width = 0.1)) + 
    xlab("Has a nature name") + ylab("Mean student reading score")

t.read <- t.test(cas[cas$hasNature, "read"], cas[!cas$hasNature, 
    "read"])
t.math <- t.test(cas[cas$hasNature, "math"], cas[!cas$hasNature, 
    "math"])

opts_knit$set(upload.fun = imgur_upload)  # upload all images to imgur.com

set.seed(1234)
library(ggplot2)
library(lattice)

```{r}

```

```{r basicconsole}
x <- 1:10
y <- round(rnorm(10, x, 1), 2)
df <- data.frame(x, y)
df
```

x <- 1:10
y <- round(rnorm(10, x, 1), 2)
df <- data.frame(x, y)
df

##     x    y
## 1   1 1.31
## 2   2 2.31
## 3   3 3.36
## 4   4 3.27
## 5   5 5.04
## 6   6 6.11
## 7   7 8.43
## 8   8 8.98
## 9   9 8.38
## 10 10 9.27

```{r simpleplot}
plot(x)
```

plot(x)

```{r multipleplots}
boxplot(1:10~rep(1:2,5))
plot(x, y)
```

boxplot(1:10 ~ rep(1:2, 5))

plot(x, y)

qplot(x, y, data = df)

xyplot(y ~ x)

```{r dotpointprint, results='asis', echo=FALSE}
cat("Here are some dot points\n\n")
cat(paste("* The value of y[", 1:3, "] is ", y[1:3], sep="", collapse="\n"))
```

```{r createtable, results='asis', echo=FALSE}
cat("x | y", "--- | ---", sep="\n")
cat(apply(df, 1, function(X) paste(X, collapse=" | ")), sep = "\n")
```

```{r echo=FALSE, comment="", echo=FALSE}
head(df)
```

  x    y
1 1 1.31
2 2 2.31
3 3 3.36
4 4 3.27
5 5 5.04
6 6 6.11

```{r smallplot, fig.width=3, fig.height=3}
plot(x)
```

plot(x)

```{r longanalysis, cache=TRUE}
for (i in 1:5000) {
    lm((i+1)~i)
}
```

```{r}
x <- 1:10
x
```

sed -i -e '1,/<body>/d' -e'/^<\/body>/,$d' -e 's/<br\/>$//' filename.html

``` {r optsknit}
opts_knit$set(upload.fun = imgur_upload) # upload all images to imgur.com
```

library(ggplot2)

set.seed(4444)
Data <- expand.grid(group=c("Apples", "Bananas", "Carrots", "Durians", 
            "Eggplants"),
            year=c("2000", "2001", "2002"),
            quality=c("Grade A", "Grade B", "Grade C", "Grade D", 
            "Grade E"))
Group.Weight <- data.frame(
    group=c("Apples", "Bananas", "Carrots", "Durians", "Eggplants"),
    group.weight=c(1,1,-1,0.5, 0))
Quality.Weight <- data.frame(
    quality=c("Grade A", "Grade B", "Grade C", "Grade D", "Grade E"),
    quality.weight = c(1,0.5,0,-0.5,-1))
Data <- merge(Data, Group.Weight)
Data <- merge(Data, Quality.Weight)
Data$score <- Data$group.weight + Data$quality.weight + 
    rnorm(nrow(Data), 0, 0.2)
Data$proportion.tasty <- exp(Data$score)/(1 + exp(Data$score))

ggplot(data=Data, 
       aes(x=factor(year), y=proportion.tasty, 
           group=group,
           shape=group,
           color=group)) + 
               geom_line() + 
               geom_point() +
               opts(title = 
               "Proportion Tasty by Year, Quality, and Group") +
               scale_x_discrete("Year") +
               scale_y_continuous("Proportion Tasty") + 
        facet_grid(.~quality )

#!/bin/bash
for file in * ; do mv "$file" `echo "$file" | sed 's/ /-/g'`; done
for file in * ; do mv "$file" `echo "$file" | sed 's/.*/\L&/g'`; done

`for file in * ; do mv $file `echo $file | sed 's/_//'` ; done`

sudo mkdir /etc/adobe
sudo su
sudo echo "OverrideGPUValidation = 1" >> /etc/adobe/mms.cfg"

# general stimuli
<text cross>
/position=(50%, 50%)
/items=("+")
/fontstyle = ("Arial", 60pt, false)
</text>

<text error.toosoon>
/position=(50%, 50%)
/items=("PRESSED TOO SOON")
/fontstyle = ("Arial", 30pt, false)
/txcolor =(255,0,0)
</text>

function! RunInquisit()
    let cmdstring = "!start inquisit " . '"%:p"'
    exe cmdstring
endfunction

function! InquisitLevel(elements)
    " elements: 0 or 1
    if getline(v:lnum) =~ '^# .*$'
        return ">1"
    endif
    if a:elements && getline(v:lnum) =~ '^<[a-zA-Z].*$'
        return ">2"
    endif
    return "=" 
endfunction                  
au BufEnter *.exp setlocal foldexpr=InquisitLevel(1)
au BufEnter *.exp setlocal foldmethod=expr

syntax case ignore
syn clear

syn keyword inquisitElementKeyword contained survey 
syn keyword inquisitElementKeyword contained surveypage 
syn keyword inquisitElementKeyword contained textbox 
...

syn match inquisitOperator contained "||"
syn match inquisitOperator contained "&&"

syn match inquisitNumber contained "\<[0-9%][0-9%]*"

syn region inquisitString start=+"+ end=+"+

hi link inquisitElementKeyword Keyword
hi link inquisitOperator Operator
hi link inquisitDelimiter Delimiter  
hi link inquisitString Character

au BufEnter *.otl setlocal tabstop = 2
au BufEnter *.otl setlocal shiftwidth=2

text1
    text2.1
    text2.2
    text2.3
    text2.4

au FileType vo_base  map <leader>o kjzO^y0zco<Esc>0DPA<Esc>zxA

<<initial_settings, echo=false>>=
options(stringsAsFactors=FALSE)
options(width=80)
library(psych) # used scoring and alpha
library(CTT) # used for spearman brown prophecy
@

<<import_data, echo=false>>=
cases <- read.delim("data/cases.tsv")
items <- read.delim("meta/items.tsv")
items$variable <- paste("item", items$item, sep="")
@

...
The example involves performing an item analysis of 
responses of \Sexpr{nrow(cases)} students 
to a set of \Sexpr{nrow(items)} multiple choice test items.
...

<<>>=
<<initial_settings>>
<<import_data>>
@

<<score_test>>=
itemstats <- score.multiple.choice(key = items$correct, 
            data = cases[,items$variable])
@

<<plot_mean_by_r, fig=true>>=
plot(r ~ mean , itemstats$item.stats, type="n")
text(itemstats$item.stats$mean, itemstats$item.stats$r, 1:50)
abline(h=.2, v=c(.5, .9))
@

<<flag_bad_items>>=
rules <- list(
        tooEasy = .95,
        tooHard = .3,
        lowR = .15)
oritemstats$item.stats$tooEasy <- 
    oritemstats$item.stats$mean > rules$tooEasy
...
@

\begin{itemize}
\item \emph{Too Easy}: mean correct $>$
\Sexpr{rules$tooEasy}.
\Sexpr{sum(oritemstats$item.stats$tooEasy)}
items were bad by this definition.
... 
\end{itemize}

The formula suggests  that in order to obtain
an alpha of \Sexpr{sbrown$targetAlpha},
\Sexpr{round(sbrown$multiple, 2)} times as many items are required.
Thus, the final scale would need around
\Sexpr{ceiling(sbrown$refinedItemCount)} items.
Assuming a similar number of good and bad items,
this would require an initial pool of around
\Sexpr{ceiling(sbrown$totalItemCount)} items.

	Year	Prop. Female
1	1924	0.07
2	1928	0.08
3	1932	0.08
4	1936	0.12
5	1948	0.18
6	1952	0.23
7	1956	0.26
8	1960	0.38
9	1964	0.37
10	1968	0.37
11	1972	0.36
12	1976	0.35
13	1980	0.34
14	1984	0.36
15	1988	0.37
16	1992	0.43
17	1994	0.43
18	1998	0.44
19	2002	0.45
20	2006	0.46

Rank	Country	Total	%
1	NOR	280	12.12%
2	USA	216	9.35%
3	URS	194	8.39%
4	AUT	185	8.01%
5	GER	158	6.84%
6	FIN	151	6.53%
7	CAN	119	5.15%
8	SUI	118	5.11%
9	SWE	118	5.11%
10	GDR	110	4.76%

Year	City	Discipline	Event	Medal
1994	Lillehammer	Short Track S.	5000m relay	Bronze
1998	Nagano	Alpine Skiing	slalom	Bronze
2002	Salt Lake City	Short Track S.	1000m	Gold
2002	Salt Lake City	Freestyle Ski.	aerials	Gold
2006	Turin	Freestyle Ski.	aerials	Bronze
2006	Turin	Freestyle Ski.	moguls	Gold

Year	Gold	Silver
1924	CAN	USA
1928	CAN	SWE
1932	CAN	USA
1936	GBR	CAN
1948	CAN	TCH
1952	CAN	USA
1956	URS	USA
1960	USA	CAN
1964	URS	SWE
1968	URS	TCH
1972	URS	USA
1976	URS	TCH
1980	USA	URS
1984	URS	TCH
1988	URS	FIN
1992	EUN	CAN
1994	SWE	CAN
1998	CZE	RUS
2002	CAN	USA
2006	SWE	FIN

A	B	C
1	Male	Blue
2	Female	Pink

Jeromy Anglim's Blog: Psychology and Statistics

Google Reader Replacements: Feedly and The Old Reader

Demise of Google Reader

Feedly

Alternatives

Beamer presentations using pandoc, markdown, LaTeX, and a makefile

Motivation

Existing resources on combining Beamer, markdown and pandoc

My approach

makefile

Preamble LaTeX file: main.tex

Body markdown slide file: talk.md

Conclusion

Other aspects

Syntax highlighting of markdown+LaTeX in Vim

Code folding of markdown-Beamer

Showing backticks and single quotes properly in code

Video: knitr, R Markdown, and R Studio: Introduction to Reproducible Analysis

Talk Overview

Video of Talk

More Videos from Melbourne R Users

Relevant links:

Converting Sweave LaTeX to knitr LaTeX: A case study

Additional Resources

My conversion from Sweave to knitr

Additional Resources

How to Convert Sweave LaTeX to knitr R Markdown: Winter Olympic Medals Example

Overview

Converting from Sweave to R Markdown

R code chunks

Figures and Tables

Basic formatting

LaTeX things

R Markdown Analysis of Winter Olympic Medal Data

Import Dataset

Total Medals by Year

Gender Ratio by Year

Countries with the Most Medals

Proportion of Gold Medals by Country

How many different countries have won medals by year?

Australia at the Winter Olympics

Ice Hockey

Reflections on the Conversion Process

Additional Resources

Example Reproducible Report using R Markdown: Analysis of California Schools Test Data

Preliminaries

Load packages and data

Q 1 What does the CASchools dataset involve?

Q. 2 To what extent does expenditure per student vary?

Q. 3a What predicts expenditure per student?

Q. 4 what is the relationship between district level maths and reading scores?

Q. 5 What is the relationship between maths and reading after partialling out other effects?

Q. 6 What fraction of a computer does each student have?

Q. 7 What is a good model of the combined effect of other variables on academic performance (i.e., math and read)?

Now, let's look at some weird stuff.

Q. 8.1 What are common words in Californian School names?

Q. 8.2 Is the number of letters in the school's name related to academic performance?

Q. 8.3 Is the number of words in the school name related to academic performance?

Q. 8.4 Are schools with nice popular nature words in their name doing better academically?

Review

Getting Started with R Markdown, knitr, and Rstudio 0.96

Getting started

Prepare for analyses

Basic console output

Plots

Simple plot

Multiple plots

ggplot2 plot

lattice plot

R Code chunk features

Create Markdown code from R

Create Markdown table code from R

Control output display

Control figure size

Cache analysis

Basic markdown functionality

Dot Points

Equations

Tables

Hyperlinks

Preamble LaTeX file: `main.tex`

Body markdown slide file: `talk.md`

`ggplot2` plot

`lattice` plot

Annoying `<br/>`'s