Saturday, July 26, 2014

Guns are Cool - States

Last week I looked at time effects of the shootingtracker database. This week I will look at the states. Some (smaller) states never made it on the database. Other states, far too frequently. The worst of these California. After correcting for population size, Washington DC with 0.8 shootings per year per 100000 inhabitants sticks out. Subsequently, the observed incidences have been shrunk towards the common distribution. In these posterior estimates Massachusetts is best while Washington DC, Illinois and Louisiana are worst.

Data

Data are from the shootingtracker database and have been updated July 25. The last entry in the database is July 25. The script is a slightly adapted version from last week. Most notably, I did not notice last week the database uses KS for Kansas, but KA is what my other data uses. Note that in order make this post readable, all code is at a separate section.

Adding state information

The first thing is to get one record per state. The datasets package has state information, while number of inhabitants (2013 data) is plucked from census.gov. After merging some states have missing counts for shootings absence in the shootingtracker database is taken as 0 shootings.
Incidence is calculated per 100 000 inhabitants per year (365 days).

Calculating posteriors

Assuming a beta prior and binomial distribution for number of shootings, Jags is used to estimate posteriors. The one noticeable thing in the model is the prior for b. b really wants to be high, hence the prior is quite a bit wider than usual. The current values seem to be acceptable, see also the output at the bottom of this post.
The plot shows a number of things. In red are the observed values. In black the posterior values; means and 90% intervals of the samples. The blue line is the value to which the data is shrunk. Washington DC is shrunk quite a bit, because it has relatively few people and a large value. The same is true for Delaware. Illinois on the other hand has many inhabitants, hence is shrunk much less. It also has a more narrow interval. On the other side, states like Wyoming are shrunk upwards. Again, more inhabitants means less shrinking, Massachusetts ends up as having least number of shootings. Really big states, New York, Texas and California, are not shrunk at all and have quite narrow intervals. Texas and NY on the good side, California on the bad side. 



The map shows the posteriors in a different way. At this point it might be worth noting that the north east is actually doing quite well, but since these states have so little inhabitants they are shrunk upwards and look bad. The current approach does not provide a natural way to avoid this. In a year or two additional data may cause at least some change in posteriors of the smaller states, both in the top and the bottom of this list.

Code

reading data

r13 <- readLines('raw13.txt')
r14 <- readLines('raw14.txt')
r1 <- c(r13,r14)
r2 <- gsub('\\[[a-zA-Z0-9]*\\]','',r1)
r3 <- gsub('^ *$','',r2)
r4 <- r3[r3!='']
r5 <- gsub('\\t$','',r4)
r6 <- gsub('\\t References$','',r5)
r7 <- read.table(textConnection(r6),
    sep='\t',
    header=TRUE,
    stringsAsFactors=FALSE)
r7$Location[r7$Location=='Washington DC'] <-
    'WashingtonDC, DC'
r8 <- read.table(textConnection(as.character(r7$Location)),
    sep=',',
    col.names=c('Location','State'),
    stringsAsFactors=FALSE)
r8$State <- gsub(' ','',r8$State)
r8$State[r8$State=='Tennessee'] <- 'TN'
r8$State[r8$State=='Ohio'] <- 'OH'
r8$State[r8$State %in% c('Kansas','KA')] <- 'KS'
r8$State[r8$State=='Louisiana'] <- 'LA'
r8$State[r8$State=='Illinois'] <- 'IL'
r8$State <- toupper(r8$State)
r7$StateAbb <- r8$State
r7$Location <- r8$Location
r7 <- r7[! (r7$State %in% c( 'PUERTORICO','PR')),]
r7$Date <- gsub('/13$','/2013',r7$Date)
r7$date <- as.Date(r7$Date,format="%m/%d/%Y")

state information

r7$one <- 1
agg1 <- with(r7,
    aggregate(one,by=list(StateAbb=StateAbb),FUN=sum))
library(datasets)
states <- data.frame(StateAbb=as.character(state.abb),
    StateArea=state.area,
    State=as.character(state.name),
    StateRegion=state.region,
    Frost=state.x77[,'Frost']
)
# http://www.census.gov/popest/data/state/totals/2013/index.html
# data pre processed so the csv can be used as is.
inhabitants <- read.csv('NST-EST2013-01.treated.csv')
#put it all together
states <- rbind(states,data.frame(StateAbb='DC',
        StateArea=68.3,
        State='District of Columbia',
       StateRegion= 'South',Frost=NA))
r9 <- merge(states,inhabitants)
r10 <- merge(r9,agg1,all=TRUE)
# No data in some states
r10$x[is.na(r10$x)] <- 0
r10$incidence <- 100000*r10$x/r10$Population*365/

        as.numeric((max(r7$date)-min(r7$date)))

Jags model


library(R2jags)
datain <- list(
    count=r10$x,
    Population = r10$Population,
    n=nrow(r10),
    scale=100000*365/
        as.numeric((max(r7$date)-min(r7$date))))

model1 <- function() {
  for (i in 1:n) {
    count[i] ~ dbin(p[i],Population[i])
    p[i] ~ dbeta(a,b)
    pp[i] <- p[i]*scale
  }
  a ~ dnorm(0,.000001) %_% T(0,)
  b ~ dnorm(1e7,.00000000000001) %_% T(0,)
  betamean <- scale * a/(a+b)
}

parameters <- c('a','b','betamean','pp')
jagsfit <- jags(datain, 
    model=model1, 
    parameters=parameters,
    n.iter=20000,
    n.burnin=1000,
    n.chain=4,
    progress.bar="gui")
jagsfit

Plot 1

samples <- jagsfit$BUGSoutput$sims.matrix[,5:55]
r10$pred <- apply(samples,2,mean)
r10$l05 <- apply(samples,2,function(x) quantile(x,.05))
r10$l95 <- apply(samples,2,function(x) quantile(x,.95))
library(ggplot2)
r11 <- subset(r10,,c(State,l05,l95,pred,incidence))
r11 <- r11[order(r11$pred),]
r11$State <- factor(r11$State,levels=r11$State)
p <- ggplot(r11, aes(y=pred, x=State))
limits <- aes(ymin =l05, ymax=l95)
dodge <- position_dodge(width=0.9)
p +geom_point(position=dodge, stat="identity") +
geom_errorbar(limits, position=dodge, width=0.25) +
theme(legend.position = "bottom") +
geom_point(aes(y=incidence,x=State),
position=dodge, stat="identity",col='red') +
ylab('Incidence') +
geom_hline(yintercept=mean(jagsfit$BUGSoutput$sims.matrix[,3]),col='blue') +
coord_flip()

Plot 2

# http://stackoverflow.com/questions/20434343/color-us-states-according-to-value-in-r
# your values (here as named vector)

library(maps)
values <- (r10$pred-min(r10$pred))/(max(r10$pred)-min(r10$pred))
#values <- (r10$incidence-min(r10$incidence))/
#    (max(r10$incidence[r10$incidence<.5])-min(r10$incidence))

names(values) <- r10$state.abb
# getting the names used by map
tmp <- map('state',plot=FALSE,namesonly=TRUE)
# matching (after adjusting using gsub and tolower)
tmp <- match(gsub('(:.*)','',tmp),tolower(r10$State))
map('state',
    fill=TRUE,
    col=rgb(colorRamp(c('green','red'))(values),max=255)[tmp])

Jags output

Inference for Bugs model at "C:/Users/Kees/AppData/Local/Temp/RtmpwbAFvR/modela3419a366f0.txt", fit using jags,
 4 chains, each with 20000 iterations (first 1000 discarded), n.thin = 19
 n.sims = 4000 iterations saved
             mu.vect     sd.vect        2.5%         25%         50%
a              9.732       6.530       2.931       5.662       7.924
b        6207140.792 4125444.443 1898003.649 3609272.959 5094169.482
betamean       0.101       0.008       0.087       0.096       0.101
pp[1]          0.087       0.032       0.032       0.064       0.085
pp[2]          0.123       0.030       0.073       0.103       0.120
pp[3]          0.071       0.026       0.026       0.053       0.070
pp[4]          0.099       0.024       0.058       0.083       0.098
pp[5]          0.130       0.014       0.105       0.120       0.130
pp[6]          0.081       0.023       0.041       0.065       0.079
pp[7]          0.097       0.027       0.050       0.078       0.094
pp[8]          0.130       0.041       0.066       0.101       0.123
pp[9]          0.110       0.017       0.079       0.097       0.109
pp[10]         0.094       0.020       0.059       0.081       0.093
pp[11]         0.079       0.030       0.027       0.058       0.077
pp[12]         0.071       0.026       0.026       0.052       0.069
pp[13]         0.076       0.029       0.024       0.055       0.074
pp[14]         0.177       0.029       0.127       0.158       0.176
pp[15]         0.122       0.027       0.077       0.103       0.120
pp[16]         0.122       0.033       0.069       0.099       0.118
pp[17]         0.109       0.028       0.062       0.089       0.106
pp[18]         0.159       0.037       0.100       0.132       0.154
pp[19]         0.056       0.020       0.022       0.041       0.054
pp[20]         0.087       0.023       0.048       0.071       0.086
pp[21]         0.090       0.031       0.035       0.068       0.087
pp[22]         0.124       0.024       0.083       0.107       0.122
pp[23]         0.080       0.023       0.040       0.063       0.078
pp[24]         0.127       0.028       0.080       0.108       0.124
pp[25]         0.079       0.026       0.032       0.061       0.078
pp[26]         0.083       0.031       0.029       0.062       0.081
pp[27]         0.121       0.023       0.080       0.105       0.119
pp[28]         0.087       0.032       0.029       0.065       0.086
pp[29]         0.092       0.030       0.041       0.072       0.090
pp[30]         0.080       0.031       0.027       0.057       0.078
pp[31]         0.106       0.022       0.068       0.090       0.104
pp[32]         0.108       0.032       0.055       0.086       0.104
pp[33]         0.139       0.037       0.081       0.113       0.135
pp[34]         0.076       0.014       0.051       0.066       0.075
pp[35]         0.097       0.019       0.063       0.084       0.096
pp[36]         0.123       0.031       0.074       0.101       0.120
pp[37]         0.071       0.024       0.029       0.054       0.070
pp[38]         0.094       0.018       0.062       0.081       0.093
pp[39]         0.106       0.035       0.049       0.082       0.102
pp[40]         0.110       0.027       0.063       0.091       0.108
pp[41]         0.086       0.032       0.030       0.064       0.084
pp[42]         0.122       0.026       0.078       0.103       0.119
pp[43]         0.078       0.013       0.055       0.069       0.077
pp[44]         0.081       0.027       0.034       0.062       0.080
pp[45]         0.111       0.023       0.070       0.094       0.109
pp[46]         0.089       0.034       0.033       0.066       0.086
pp[47]         0.081       0.021       0.043       0.066       0.080
pp[48]         0.065       0.021       0.029       0.050       0.064
pp[49]         0.103       0.032       0.051       0.080       0.100
pp[50]         0.090       0.034       0.032       0.066       0.087
pp[51]         0.185       0.063       0.095       0.139       0.174
deviance     242.511      13.874     217.150     232.666     241.477
                 75%        97.5%  Rhat n.eff
a             11.529       28.847 1.068    63
b        7337119.577 18410168.237 1.067    64
betamean       0.106        0.117 1.001  4000
pp[1]          0.106        0.158 1.003   910
pp[2]          0.140        0.190 1.001  2600
pp[3]          0.088        0.125 1.006   460
pp[4]          0.114        0.152 1.001  4000
pp[5]          0.139        0.159 1.002  1800
pp[6]          0.096        0.130 1.005   560
pp[7]          0.112        0.157 1.001  4000
pp[8]          0.151        0.230 1.001  3000
pp[9]          0.120        0.147 1.001  2800
pp[10]         0.107        0.137 1.003  1100
pp[11]         0.096        0.143 1.007   450
pp[12]         0.087        0.126 1.005   550
pp[13]         0.094        0.137 1.007   420
pp[14]         0.196        0.238 1.012   250
pp[15]         0.138        0.182 1.001  4000
pp[16]         0.140        0.200 1.001  2900
pp[17]         0.124        0.171 1.001  4000
pp[18]         0.182        0.242 1.008   350
pp[19]         0.069        0.098 1.008   320
pp[20]         0.102        0.137 1.002  2000
pp[21]         0.108        0.158 1.003  1200
pp[22]         0.139        0.178 1.003   970
pp[23]         0.094        0.128 1.002  1800
pp[24]         0.144        0.190 1.004   660
pp[25]         0.096        0.135 1.006   550
pp[26]         0.102        0.151 1.003  3700
pp[27]         0.135        0.170 1.002  1800
pp[28]         0.106        0.159 1.004  1200
pp[29]         0.110        0.160 1.004  1200
pp[30]         0.099        0.145 1.005   520
pp[31]         0.120        0.154 1.002  1400
pp[32]         0.127        0.179 1.003  1800
pp[33]         0.161        0.224 1.006   450
pp[34]         0.085        0.105 1.006   480
pp[35]         0.108        0.138 1.001  3600
pp[36]         0.141        0.191 1.002  1300
pp[37]         0.087        0.123 1.004   700
pp[38]         0.106        0.133 1.001  4000
pp[39]         0.125        0.184 1.002  4000
pp[40]         0.126        0.171 1.002  4000
pp[41]         0.106        0.154 1.003  1200
pp[42]         0.138        0.180 1.006   470
pp[43]         0.086        0.106 1.006   480
pp[44]         0.098        0.143 1.005   590
pp[45]         0.125        0.161 1.002  2500
pp[46]         0.109        0.164 1.004   700
pp[47]         0.094        0.126 1.004   700
pp[48]         0.080        0.110 1.008   320
pp[49]         0.120        0.173 1.004   790
pp[50]         0.109        0.166 1.003  1400
pp[51]         0.217        0.337 1.012   210
deviance     251.840      271.487 1.039    88

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 92.8 and DIC = 335.3
DIC is an estimate of expected predictive error (lower deviance is better).

Saturday, July 19, 2014

Guns are cool - time effects

September last year I made a post using the shootingtracker data. It is attempted in shootingtracker to register all shootings with at least four victims, be they wounded or dead. The data starts  January 1st 2013, which means that by now the amount of data has almost doubled. This surely is a dataset where I hope the makers find less and less data to add. Analysis shows Sundays in summer have the highest number of shootings on a day. Three to four shootings on Sunday in July and August.

Data

Data sits in two pages on shootingtracker, 2013 and 2014. In preparation for this post I copied/pasted those data in notepad and removed headers and footers. In the 2013 data I kept names of columns. The first steps of reading the data are removing some of the (for me) extraneous info, such as reference where the data came from. Subsequently the state and town are separated and a few records which do not have the correct state abbreviation are corrected. Finally, 13 is reformatted to 2013 and date is created. The last record used is from July 9th, 2014.
r13 <- readLines('raw13.txt')
r14 <- readLines('raw14.txt')
r1 <- c(r13,r14)
head(r1)
[1] "Number\t Date\t Alleged Shooter\t Killed\t Wounded\t Location\t References"
[2] "1\t1/1/13\tCarlito Montoya\t4\t0\tSacramento, CA\t"                        
[3] " [Expand] "                                                          
[4] "2\t1/1/13\tUnknown\t1\t3\tHawthorne, CA\t"                                 
[5] " [Expand] "                                                          
[6] "3\t1/1/13\tJulian Sims\t0\t4\tMcKeesport, PA\t"
tail(r1)
[1] "141\t7/8/2014\tUnknown\t1\t4\tSan Bernardino, CA\t"    
[2] " [Expand] "                                      
[3] "142\t7/8/2014\tUnknown\t0\t5\tProvidence, RI\t"        
[4] " [Expand] "                                      
[5] "143\t7/9/2014\tRonald Lee Haskell\t6\t1\tHouston, TX\t"
[6] " [Expand] " 
r2 <- gsub('\\[[a-zA-Z0-9]*\\]','',r1)
r3 <- gsub('^ *$','',r2)
r4 <- r3[r3!='']
r5 <- gsub('\\t$','',r4)
r6 <- gsub('\\t References$','',r5)
r7 <- read.table(textConnection(r6),
    sep='\t',
    header=TRUE,
    stringsAsFactors=FALSE)
r7$Location[r7$Location=='Washington DC'] <-
    'Washington, DC'
r8 <- read.table(textConnection(as.character(r7$Location)),
    sep=',',
    col.names=c('Location','State'),
    stringsAsFactors=FALSE)
r8$State <- gsub(' ','',r8$State)
r8$State[r8$State=='Tennessee'] <- 'TN'
r8$State[r8$State=='Ohio'] <- 'OH'
r8$State[r8$State=='Kansas'] <- 'KS'
r8$State[r8$State=='Louisiana'] <- 'LA'
r8$State[r8$State=='Illinois'] <- 'IL'
r8$State <- toupper(r8$State)

r7$State <- r8$State
r7$Location <- r8$Location
r7 <- r7[r7$State != 'PUERTORICO',]
Sys.setlocale(category = "LC_TIME", locale = "C")
r7$Date <- gsub('/13$','/2013',r7$Date)
r7$date <- as.Date(r7$Date,format="%m/%d/%Y")

Effect of day and month

Effect of day of the week is pretty easy to plot, just add the day and run qplot(weekday). In this case complexities arise because I want the days in a specific order, Monday to Sunday. Second is that not all days occur equally often in the test period. This is not enough to invalidate the plot, but since I had to correct for occurrence of months for a similar plot, I decided to reuse that code for weekdays. The data.frame alldays is used to calculate the number of days in the data set. I am not going to over analyze this, Sundays stick out in a negative way.
library(ggplot2)
r7$weekday <- factor(format(r7$date,'%a'),
    levels=format(as.Date('07/07/2014',format="%m/%d/%Y")+0:6,'%a'))
r7$month <- factor(format(r7$date,'%b'),
    levels=format(as.Date('01/15/2014',format="%m/%d/%Y")+
            seq(0,length.out=12,by=30),'%b'))
alldays <- data.frame(
    date=seq(min(r7$date),max(r7$date),by=1))
alldays$weekday <- factor(format(alldays$date,'%a'),
    levels=levels(r7$weekday))
alldays$month <- factor(format(alldays$date,'%b'),
    levels=format(as.Date('01/15/2014',format="%m/%d/%Y")+
            seq(0,length.out=12,by=30),'%b'))
ggplot(data=data.frame(prop=as.numeric(table(r7$weekday)/
        table(alldays$weekday)),
    weekday=factor(levels(r7$weekday),
        levels=levels(r7$weekday))),
    aes(y=prop,x=weekday)) +
    geom_bar(stat='identity') +
    ylab('Shootings per day') +
    xlab('Day of the week')
In terms of months, it seems summer is worse than the other seasons and winter is best.
ggplot(data=data.frame(prop=as.numeric(table(r7$month)/
                    table(alldays$month)),
            month=factor(levels(r7$month),
                levels=levels(r7$month))),
        aes(y=prop,x=month)) +
    geom_bar(stat='identity') +
    ylab('Shootings per day') +
    xlab('Month')   

A model

It is not obvious, given the unequal distributions of weekdays over months, how significant a month effect is. To examine this, I have reorganized the data to display shootings per day. Data frame alldays is used again, now to ensure data with no shootings are correctly represented. The modeling shows a clear effect of months and the interaction of days and months on the brink of significance.
r7$one <- 1
ag <- aggregate(r7$one,
    by=list(date=r7$date),FUN=sum)
counts <- merge(alldays,ag,all=TRUE)
counts$x[is.na(counts$x)] <- 0
g0 <- glm(x ~ weekday ,data=counts,
    family='poisson')
g1 <- glm(x ~ weekday + month,data=counts,
    family='poisson')
g2 <- glm(x ~ weekday * month,data=counts,
    family='poisson')
anova(g0,g1,g2,test='Chisq')
Analysis of Deviance Table

Model 1: x ~ weekday
Model 2: x ~ weekday + month
Model 3: x ~ weekday * month
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
1       548     653.86                        
2       537     624.86 11   29.008 0.002264 **
3       471     538.99 66   85.870 0.050718 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Predictions

To understand the interaction, a plot is made of the expected values by day and month. Looking at this plot, the Sunday effect is most pronounced in Summer. In addition I would not be surprised if next Sunday has four shootings again.
pred1 <- expand.grid(weekday=factor(levels(r7$weekday),
        levels=levels(r7$weekday)),
    month=factor(levels(r7$month),
        levels=levels(r7$month)))
preds <- predict(g2,pred1,
    type='response',
    se=TRUE)
pred1$fit <- preds$fit
pred1$se.fit <- preds$se.fit

limits <- aes(ymax = fit+se.fit, ymin=fit-se.fit)
p <- ggplot(pred1, aes(fill=weekday, y=fit, x=month))
dodge <- position_dodge(width=0.9)
p + geom_bar(position=dodge, stat="identity") +
    geom_errorbar(limits, position=dodge, width=0.25) +
    theme(legend.position = "bottom") +
    ylab('Shootings per Day') +
    guides (fill=guide_legend('Day'))


Follow the example?

One might ask if shootings which have much media attention give cause to copycats. This is not easy to analyze, given that clearly time effects from day and month exist. Besides, which shootings get a lot of media attention? Yet we can look at the number of shootings over time and at least add the shootings with most victims in the plot. The number of victims to be marked has arbitrarily been chosen as 18 or more. In this plot I cannot see the connection. 
r7$Victims <- r7$Killed+r7$Wounded
table(r7$Victims)
  4   5   6   7   8   9  10  12  13  14  18  19  20  21 
322 104  25  26  11   5   4   2   2   2   1   1   1   1 

p <- ggplot(counts, aes(y=x,x=date))
p + stat_smooth(span=.11) +
    geom_point() +
    geom_vline(xintercept=as.numeric(r7$date[r7$Victims>15])) +
    ylab('shootings per day')

An ACF

Mostly because making an ACF is integral part of analyzing time series. Here it is, in case anybody doubted a week effect. I chose this fairly long lag, because 10 weeks seemed nice to me. 
plot(acf(counts$x,lag.max=70))

Saturday, July 12, 2014

odfweave setup and counting logicals

Two short items in this blogpost. Since it was not obvious how to run odfWeave() in my particular setup, the call I am using. Then there were several people crosstabulating logical vectors, so I wanted to play along, 80 times faster than table().

odfWeave

My particular setup consists of R, 7-zip, libreoffice. Somehow they don't 100% play along when using odfWeave. I had that problem this spring and decided to put my solution in a post at some point. In terms of versions therefore, I had that with my previous versions, and tested that it still runs with my new setup (R 3.1.1, Libreoffice  4.2.5.2). The only loose end, is that odfWeave complains I am re-using a directory, and that I need to empty said directory manually.
# the standard example call that works for me
demoFile <- system.file("examples", "simple.odt", package = "odfWeave")
outputFile <- gsub("simple.odt", "output.odt", demoFile)
odfWeave(demoFile, outputFile,
    workDir='C:\\Users\\Kees\\Documents\\tmp',
    odfWeaveControl(zipCmd = 
            c("C:\\Progra~1\\7-Zip\\7z a -tzip $$file$$ . -r", 
                "C:\\Progra~1\\7-Zip\\7z x -tzip $$file$$ -yr") ))
# removing files
file.remove(dir('C:\\Users\\Kees\\Documents\\tmp',
        recursive=TRUE,
        full.names=TRUE))

# using a different directory
odfWeave('C:\\Users\\Kees\\Documents\\test\\testcases.odt',
    'C:\\Users\\Kees\\Documents\\test\\testout.odt',
    workDir='C:\\Users\\Kees\\Documents\\tmp',
    odfWeaveControl(zipCmd = 
            c("C:\\Progra~1\\7-Zip\\7z a -tzip $$file$$ . -r", 
                "C:\\Progra~1\\7-Zip\\7z x -tzip $$file$$ -yr") ))

Cross table of logical vectors

This was started in Sometimes Table is not the Answer – a Faster 2×2 Table and carried on with Sometimes I feel (some) need for speed. So, I wanted to add my own attempts. The aim is to make a cross table of two logical vectors with a minimum of time. Which becomes important if these vectors are long. Solutions from previous posts.
set.seed(2014)

manual = sample(c(TRUE, FALSE), 10e6, replace = TRUE)
auto = sample(c(TRUE, FALSE), 10e6, replace = TRUE)

logical.tab = function(x, y) {
  tt = sum(x & y)
  tf = sum(x & !y)
  ft = sum(!x & y)
  ff = sum(!x & !y)
  return(matrix(c(ff, tf, ft, tt), 2, 2))
}

basic.tab2 = function(x, y) {
  dif = x - y
  tf = sum(dif > 0)
  ft = sum(dif < 0)
  tt = sum(x*y)
  ff = length(dif) - tt - tf - ft
  return(c(tf, ft, tt, ff))
}
tabulate(manual + auto *2+1, 4)

My idea was we should use the margins and go back from there.
my.tab = function(x, y) {
  tt = sum(x * y)
  t1=sum(x)
  t2=sum(y)
  return(matrix(c(length(x)-t1-t2+tt,  t1-tt, t2-tt, tt), 2, 2))
}

my.tab2 <- function(x, y) {
  phase1 <- colSums(cbind(x,y,x*y))
  return(matrix(c(length(x)-sum(phase1[-3])+phase1[3],
     phase1[-3]-phase1[3],
     phase1[3]),2,2))
}
With my particular hardware table() is just too slow to microbenchmark often, but 80 times faster than table() is not bad.
library(microbenchmark)
microbenchmark(
    logical.tab(manual, auto), 
    basic.tab2(manual, auto),
    my.tab(manual,auto),
    my.tab2(manual,auto),
    tabulate(manual + auto *2+1, 4),
    table(manual,auto),
    times = 20)
Unit: milliseconds
                               expr        min         lq     median         uq        max neval
          logical.tab(manual, auto)  2852.5587  2888.8590  2906.4571  2972.3916  3227.0821    20
           basic.tab2(manual, auto)   705.8153   722.5800   746.1683   765.9400   957.5435    20
               my.tab(manual, auto)   185.8359   186.6829   188.0988   224.2308   413.5623    20
              my.tab2(manual, auto)   463.2731   481.8843   487.7825   512.2563   694.1729    20
 tabulate(manual + auto * 2 + 1, 4)   276.1837   300.8009   315.9451   379.7302   534.7997    20
                table(manual, auto) 15703.0576 16132.0100 16231.3342 16466.7445 19012.0273    20

Sunday, July 6, 2014

Stone Flakes V, networks again

Last week I tried pcalg. This week deal (Learning Bayesian Networks with Mixed Variables). The aim n this post I want to try something new, a causal graphical model. The aim here is just as much to get myself a feel what these things do as to understand how the stone flakes data fit together.

Data

Data are stone flakes data which I analyzed previously. The first post was clustering, second linking to hominid type, third regression. Together these made for the bulk of a standard analysis. In this new analysis the same starting data is used.
r2 <- read.table('StoneFlakes.txt',header=TRUE,na.strings='?')
r1 <- read.table('annotation.txt',header=TRUE,na.strings='?')
r12 <- merge(r1,r2) 

r12$group <- factor(r12$group,labels=c('Lower Paleolithic',
        'Levallois technique',
        'Middle Paleolithic',
        'Homo Sapiens'))
r12$site <- factor(c('other','gravel pit')[r12$site+1])
r12$mat <- factor(c('flint','other')[r12$mat])
r12$lmage <- log10(-r12$age)

Deal

The starting point of this post was to continue/ my analysis of last week. But when I discovered deal could be used to discover the model, repeating last week's analysis was chosen instead. Deal does not have a Vignette, but there is a paper
deal: A Package for Learning Bayesian Networks which helped me very well to get started.

First Model

Initially I wanted to start with a model containing only continuous variables, similar to before, but that threw an error in jointprior(). Hence I added groups as factor. Autosearch() and heuristicsearch() give a lot of output, basically one line for each step. For brievety these are not shown. The good thing about this model is that it has a solution where 'group' is driving other variables.
library(deal)
rfin <- subset(r12,,c(names(r2)[-1],'group'))
rfin <- rfin[complete.cases(rfin),]
rfin.nw <- network(rfin)
rfin.prior <- jointprior(rfin.nw)

Imaginary sample size: 8 
rfin.nw <- learn(rfin.nw,rfin,rfin.prior)$nw
rfin.search <- autosearch(rfin.nw,
    rfin,
    rfin.prior,
    trace=FALSE)

plot(rfin.search$nw)
Heuristic is used to further improve the model. In the end the model seems a bit more complex that pcalg, but not unreasonably so.
rfin.heuristic <- heuristic(rfin.search$nw,
    rfin,
    rfin.prior,
    restart=10,
    trace=FALSE,
    trylist=rfin.search$trylist)
plot(rfin.heuristic$nw)

Second model

For brevity I won't be repeating the code. It is all the same except for the data going in, which will be shown. The second model is similar to the first, but the (potential) outliers have been removed. This model looks even more clean.
rfin <- subset(r12,
    !(r12$ID %in% c('ms','c','roe','sz','va','arn')),
    c(names(r2)[-1],'group'))
rfin <- rfin[complete.cases(rfin),]

Third model

Including all sensible factors is the model I wanted to do. However, it seemed the imaginary sample size grew to 96. It is my experience that higher imaginary sample sizes produce more complex networks and longer run times. The current model seems a bit too complex to my liking. Moving under the recommended number gave runtime errors. 

Final model

In the final model (a few are skipped now) a number of simplifications were made. Region is removed as factor, log(-age) as continuous variable. Group has lost Homo Sapiens, since that category had only three records. The model is restricted in the sense that group cannot be the result of other variables. In the plot these are shown as red arrows.
rfin <- subset(r12,
    !(r12$ID %in% c('ms','c','roe','sz','va','arn')),
    c(-ID,-number,-age,-dating,-region,-lmage))
rfin <- rfin[rfin$group !='Homo Sapiens',]
rfin <- rfin[complete.cases(rfin),]
rfin.nw <- network(rfin)
rfin.prior <- jointprior(rfin.nw)
mybanlist <- matrix(
    c(2:11,
        rep(1,10)),ncol=2)
banlist(rfin.nw) <- mybanlist

Conclusion

Deal makes too complex networks for my liking, pcalg cannot use discrete variables. Deal has banned links, a feature which helps. pcalg made more nice plots, but I have the feeling that is relatively easily remedied. Neither gave a model which struck me as a model to continue with. I'll be needing quite some more study to feel comfortable with this kind of models.

Sunday, June 29, 2014

stone flakes IV

In this post I want to try something new, a causal graphical model. The aim here is just as much to get myself a feel what these things do as to understand how the stone flakes data fit together.

Data

Data are stone flakes data which I analyzed previously. The first post was clustering, second linking to hominid type, third regression. Together these made for the bulk of a standard analysis. In this new analysis the same starting data is used.
r2 <- read.table('StoneFlakes.txt',header=TRUE,na.strings='?')
r1 <- read.table('annotation.txt',header=TRUE,na.strings='?')
r12 <- merge(r1,r2)

Packages

The main package used is pcalg (Methods for graphical models and causal inference). Even though it lives on cran, it requires RBGL (An interface to the BOOST graph library) which lives on Bioconductor. Plots are made via Rgraphvis (Provides plotting capabilities for R graph objects), Bioconductor again, which itself has the hard work done by graphviz, which, on my linux machine, is a few clicks to install. 
library('pcalg')
library('Rgraphviz')


First analysis

I am just following the vignette here, to get some working code.
rx <- subset(r12,,names(r2)[-1])
rx <- rx[complete.cases(rx) & !(r12$ID %in% c('ms','c','roe','sz','va','arn')),]
suffStat <- list(C = cor(rx), n = nrow(rx))
pc.gmG <- pc(suffStat, indepTest = gaussCItest,
    p = ncol(rx), alpha = 0.01)
png('graph1.png')
plot(pc.gmG, main = "")


personally I dislike this plot since you have to know which variable is which number. I don't think that is acceptable for things one wants to share. Since I could not find documentation how to modify this via the plot statement, I took the ugly road of directly modifying an S4 object; pc.Gmc.
pc.gmG@graph@nodes <- names(rx)
names(pc.gmG@graph@edgeL) <- names(rx)
png('graph2.png')
plot(pc.gmG, main = "")
dev.off()

This makes some sense looking at the variable names.
RTI (Relative-thickness index of the striking platform) is connected to WDI (Width-depth index of the striking platform). PSF (platform primery (yes/no, relative frequency)) is related to FSF (Platform facetted (yes/no, relative frequency)). PSF is also related to PROZD (Proportion of worked dorsal surface (continuous)) which then goes to ZDF1 (Dorsal surface totally worked (yes/no, relative frequency)). ZDF1 is also influenced by FLA (Flaking angle (the angle between the striking platform and the splitting surface)).

Second analysis

Much as like this analysis, it does not lead to a connection between flakes on one hand and age or group on the other hand. Since the algorithm assumes normal distributed variables, group is out of the question. Log(-age) seems to be closest to normal distributed.
rx <- subset(r12,,names(r2)[-1])
rx$lmage <- log(-r12$age)
rx <- rx[complete.cases(rx) & !(r12$ID %in% c('ms','c','roe','sz','va','arn')),]
suffStat <- list(C = cor(rx), n = nrow(rx))
pc.gmG <- pc(suffStat, indepTest = gaussCItest,
    p = ncol(rx), alpha = 0.01)
pc.gmG@graph@nodes <- names(rx)
names(pc.gmG@graph@edgeL) <- names(rx)
plot(pc.gmG, main = "")


T

Adding age links the two parts, while keeping most of the previous graph unchanged. The causal link however, seems reversed, does age cause change in flakes or do changes in flakes cause age? Nevertheless, it does show a different picture than before. In linear regression FSF and LBI contributed, but there I had not removed the outliers. In this approach FSF features, but is in its turn driven by PSF. The other direct influence is ZDF1, which is now also driven by WDI.

Third analysis

It is probably pushing the limits of what is normal distributed, but there are two binairy variables, stone material (1=flint, 2=other) and site (1=gravel pit, 0=other).
rx <- subset(r12,,c('site','mat',names(r2)[-1]))
rx$lmage <- log(-r12$age)
rx <- rx[complete.cases(rx) & !(r12$ID %in% c('ms','c','roe','sz','va','arn')),]
suffStat <- list(C = cor(rx), n = nrow(rx))
pc.gmG <- pc(suffStat, indepTest = gaussCItest,
    p = ncol(rx), alpha = 0.01)
pc.gmG@graph@nodes <- names(rx)
names(pc.gmG@graph@edgeL) <- names(rx)
plot(pc.gmG, main = "")

This makes a link from material to RTI (relative thickness) and a connection site to FLA (flaking angle), see below.
# site (1=gravel pit, 0=other)
boxplot(FLA ~ c('other','gravel pit')[site+1],
    data=r12,
    ylab='FLA',
    xlab='site')

#stone material (1=flint, 2=other)
boxplot(RTI ~ c('flint','other')[mat],
    data=r12,
    ylab='RTI',
    xlab='mat')