## Loading required package: knitr

## The basic files and libraries needed for most presentations
# creates the libraries and common-functions sections
read_chunk("../common/utility_functions.R")

require(ggplot2) #for plots
require(lattice) # nicer scatter plots
require(plyr) # for processing data.frames
require(grid) # contains the arrow function
require(biOps) # for basic image processing
require(doMC) # for parallel code
require(png) # for reading png images
require(gridExtra)
require(reshape2) # for the melt function
## To install EBImage
# source("http://bioconductor.org/biocLite.R")
# biocLite("EBImage")
require(EBImage) # for more image processing
used.libraries<-c("ggplot2","lattice","plyr","reshape2","grid","gridExtra","biOps","png","EBImage")

# start parallel environment
registerDoMC()
# functions for converting images back and forth
im.to.df<-function(in.img,out.col="val") {
  out.im<-expand.grid(x=1:nrow(in.img),y=1:ncol(in.img))
  out.im[,out.col]<-as.vector(in.img)
  out.im
}
df.to.im<-function(in.df,val.col="val",inv=F) {
  in.vals<-in.df[[val.col]]
  if(class(in.vals[1])=="logical") in.vals<-as.integer(in.vals*255)
  if(inv) in.vals<-255-in.vals
  out.mat<-matrix(in.vals,nrow=length(unique(in.df$x)),byrow=F)
  attr(out.mat,"type")<-"grey"
  out.mat
}
ddply.cutcols<-function(...,cols=1) {
  # run standard ddply command 
  cur.table<-ddply(...)
  cutlabel.fixer<-function(oVal) {
    sapply(oVal,function(x) {
      cnv<-as.character(x)
      mean(as.numeric(strsplit(substr(cnv,2,nchar(cnv)-1),",")[[1]]))
    })
  }
  cutname.fixer<-function(c.str) {
    s.str<-strsplit(c.str,"(",fixed=T)[[1]]
    t.str<-strsplit(paste(s.str[c(2:length(s.str))],collapse="("),",")[[1]]
    paste(t.str[c(1:length(t.str)-1)],collapse=",")
  }
  for(i in c(1:cols)) {
    cur.table[,i]<-cutlabel.fixer(cur.table[,i])
    names(cur.table)[i]<-cutname.fixer(names(cur.table)[i])
  }
  cur.table
}

show.pngs.as.grid<-function(file.list,title.fun,zoom=1) {
  preparePng<-function(x) rasterGrob(readPNG(x,native=T,info=T),width=unit(zoom,"npc"),interp=F)
  labelPng<-function(x,title="junk") (qplot(1:300, 1:300, geom="blank",xlab=NULL,ylab=NULL,asp=1)+
                                        annotation_custom(preparePng(x))+
                                        labs(title=title)+theme_bw(24)+
                                        theme(axis.text.x = element_blank(),
                                              axis.text.y = element_blank()))
  imgList<-llply(file.list,function(x) labelPng(x,title.fun(x)) )
  do.call(grid.arrange,imgList)
}
## Standard image processing tools which I use for visualizing the examples in the script
commean.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    weight.sum<-sum(c.cell$weight)
    data.frame(xv=mean(c.cell$x),
               yv=mean(c.cell$y),
               xm=with(c.cell,sum(x*weight)/weight.sum),
               ym=with(c.cell,sum(y*weight)/weight.sum)
    )
  })
}

colMeans.df<-function(x,...) as.data.frame(t(colMeans(x,...)))

pca.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    c.cell.cov<-cov(c.cell[,c("x","y")])
    c.cell.eigen<-eigen(c.cell.cov)
    
    c.cell.mean<-colMeans.df(c.cell[,c("x","y")])
    out.df<-cbind(c.cell.mean,
                  data.frame(vx=c.cell.eigen$vectors[1,],
                             vy=c.cell.eigen$vectors[2,],
                             vw=sqrt(c.cell.eigen$values),
                             th.off=atan2(c.cell.eigen$vectors[2,],c.cell.eigen$vectors[1,]))
    )
  })
}
vec.to.ellipse<-function(pca.df) {
  ddply(pca.df,.(val),function(cur.pca) {
    # assume there are two vectors now
    create.ellipse.points(x.off=cur.pca[1,"x"],y.off=cur.pca[1,"y"],
                          b=sqrt(5)*cur.pca[1,"vw"],a=sqrt(5)*cur.pca[2,"vw"],
                          th.off=pi/2-atan2(cur.pca[1,"vy"],cur.pca[1,"vx"]),
                          x.cent=cur.pca[1,"x"],y.cent=cur.pca[1,"y"])
  })
}

# test function for ellipse generation
# ggplot(ldply(seq(-pi,pi,length.out=100),function(th) create.ellipse.points(a=1,b=2,th.off=th,th.val=th)),aes(x=x,y=y))+geom_path()+facet_wrap(~th.val)+coord_equal()
create.ellipse.points<-function(x.off=0,y.off=0,a=1,b=NULL,th.off=0,th.max=2*pi,pts=36,...) {
  if (is.null(b)) b<-a
  th<-seq(0,th.max,length.out=pts)
  data.frame(x=a*cos(th.off)*cos(th)+b*sin(th.off)*sin(th)+x.off,
             y=-1*a*sin(th.off)*cos(th)+b*cos(th.off)*sin(th)+y.off,
             id=as.factor(paste(x.off,y.off,a,b,th.off,pts,sep=":")),...)
}
deform.ellipse.draw<-function(c.box) {
  create.ellipse.points(x.off=c.box$x[1],
                        y.off=c.box$y[1],
                        a=c.box$a[1],
                        b=c.box$b[1],
                        th.off=c.box$th[1],
                        col=c.box$col[1])                    
}
bbox.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    c.cell.mean<-colMeans.df(c.cell[,c("x","y")])
    xmn<-emin(c.cell$x)
    xmx<-emax(c.cell$x)
    ymn<-emin(c.cell$y)
    ymx<-emax(c.cell$y)
    out.df<-cbind(c.cell.mean,
                  data.frame(xi=c(xmn,xmn,xmx,xmx,xmn),
                             yi=c(ymn,ymx,ymx,ymn,ymn),
                             xw=xmx-xmn,
                             yw=ymx-ymn
                  ))
  })
}

# since the edge of the pixel is 0.5 away from the middle of the pixel
emin<-function(...) min(...)-0.5
emax<-function(...) max(...)+0.5
extents.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    c.cell.mean<-colMeans.df(c.cell[,c("x","y")])
    out.df<-cbind(c.cell.mean,data.frame(xmin=c(c.cell.mean$x,emin(c.cell$x)),
                                         xmax=c(c.cell.mean$x,emax(c.cell$x)),
                                         ymin=c(emin(c.cell$y),c.cell.mean$y),
                                         ymax=c(emax(c.cell$y),c.cell.mean$y)))
  })
}

common.image.path<-"../common"
qbi.file<-function(file.name) file.path(common.image.path,"figures",file.name)
qbi.data<-function(file.name) file.path(common.image.path,"data",file.name)

th_fillmap.fn<-function(max.val) scale_fill_gradientn(colours=rainbow(10),limits=c(0,max.val))

Quantitative Big Imaging

author: Kevin Mader date: 23 April 2015 width: 1440 height: 900 css: ../common/template.css transition: rotate

ETHZ: 227-0966-00L

Statistics and Reproducibility

Course Outline

source('../common/schedule.R')

• 19th February - Introduction and Workflows
• 26th February - Image Enhancement (A. Kaestner)
• 5th March - Basic Segmentation, Discrete Binary Structures
• 12th March - Advanced Segmentation
• 19th March - Applying Graphical Models and Machine Learning (A. Lucchi)
• 26th March - Analyzing Single Objects
• 2nd April - Analyzing Complex Objects
• 16th April - Groups and Spatial Distribution
• 23rd April - Statistics and Reproducibility
• 30th April - Dynamic Experiments
• 7th May - Scaling Up / Big Data
• 21th May - Guest Lecture, Applications in Material Science and Wood (? and A. Patera)
• 28th May - Project Presentations

Literature / Useful References

Books

• Jean Claude, Morphometry with R
  • Online through ETHZ
  • Chapter 3
  • Buy it
• John C. Russ, “The Image Processing Handbook”,(Boca Raton, CRC Press)
  • Available online within domain ethz.ch (or proxy.ethz.ch / public VPN)
• Hypothesis Testing Chapter

Papers / Sites

• Google/Stanford Statistics Intro
  • https://www.youtube.com/watch?v=YFC2KUmEebc
• MCB 140 P-value lecture at UC Berkeley (Audio)
  • https://itunes.apple.com/us/itunes-u/mcb-140-fall-2007-general/id461120088?mt=10
• Correlation and Causation (Video)
  • https://www.youtube.com/watch?v=YFC2KUmEebc
• Matlab Unit Testing Documentation
• Visualizing Genomic Data (General Visualization Techniques)

• NIMRod Parameter Studies

• M.E. Wolak, D.J. Fairbairn, Y.R. Paulsen (2012) Guidelines for Estimating Repeatability. Methods in Ecology and Evolution 3(1):129-137.

Previously on QBI ...

• Image Enhancment
  • Highlighting the contrast of interest in images
  • Minimizing Noise
• Understanding image histograms
• Automatic Methods
• Component Labeling
• Single Shape Analysis
• Complicated Shapes
• Distribution Analysis

Quantitative "Big" Imaging

The course has covered imaging enough and there have been a few quantitative metrics, but "big" has not really entered.

What does big mean?

• Not just / even large
• it means being ready for big data
• volume, velocity, variety (3 V's)
• scalable, fast, easy to customize

So what is "big" imaging

• doing analyses in a disciplined manner
  • fixed steps
  • easy to regenerate results
  • no magic
• having everything automated
  • 100 samples is as easy as 1 sample
• being able to adapt and reuse analyses
  • one really well working script and modify parameters
  • different types of cells
  • different regions

Objectives

1. Scientific Studies all try to get to a single number
   • Make sure this number is describing the structure well (what we have covered before)
   • Making sure the number is meaningful (today!)
2. How do we compare the number from different samples and groups?
   • Within a sample or same type of samples
   • Between samples
3. How do we compare different processing steps like filter choice, minimum volume, resolution, etc?
4. How do we evaluate our parameter selection?
5. How can we ensure our techniques do what they are supposed to do?
6. How can we visualize so much data? Are there rules?

Outline

• Motivation (Why and How?)
• Scientific Goals
• Reproducibility
• Statistical metrics and results
• Parameterization
  • Parameter sweep
  • Sensitivity analysis
• Unit Testing
• Visualization

What do we start with?

Going back to our original cell image

1. We have been able to get rid of the noise in the image and find all the cells (lecture 2-4)
2. We have analyzed the shape of the cells using the shape tensor (lecture 5)
3. We even separated cells joined together using Watershed (lecture 6)
4. We have created even more metrics characterizing the distribution (lecture 7)

We have at least a few samples (or different regions), large number of metrics and an almost as large number of parameters to tune

How do we do something meaningful with it?

Correlation and Causation

One of the most repeated criticisms of scientific work is that correlation and causation are confused.

1. Correlation
   • means a statistical relationship
   • very easy to show (single calculation)
2. Causation
   • implies there is a mechanism between A and B
   • very difficult to show (impossible to prove)

Controlled and Observational

There are two broad classes of data and scientific studies.

Observational

• Exploring large datasets looking for trends
• Population is random
• Not always hypothesis driven
• Rarely leads to causation

We examined 100 people and the ones with blue eyes were on average 10cm taller

In 100 cake samples, we found a 0.9 correlation between cooking time and bubble size

Controlled

• Most scientific studies fall into this category
• Specifics of the groups are controlled
• Can lead to causation

We examined 50 mice with gene XYZ off and 50 gene XYZ on and as the foot size increased by 10%

We increased the temperature and the number of pores in the metal increased by 10%

Simple Model: Magic / Weighted Coin

incremental: true

Since most of the experiments in science are usually specific, noisy, and often very complicated and are not usually good teaching examples

• Magic / Biased Coin
  • You buy a magic coin at a shop
  • How many times do you need to flip it to prove it is not fair?
  • If I flip it 10 times and another person flips it 10 times, is that the same as 20 flips?
  • If I flip it 10 times and then multiple the results by 10 is that the same as 100 flips?
  • If I buy 10 coins and want to know which ones are fair what do I do?

Simple Model: Magic / Weighted Coin

1. Each coin represents a stochastic variable $\mathcal{X}$ and each flip represents an observation $\mathcal{X}_i$.
2. The act of performing a coin flip $\mathcal{F}$ is an observation $\mathcal{X}_i = \mathcal{F}(\mathcal{X})$

We normally assume

1. A fair coin has an expected value of $E(\mathcal{X})=0.5 \rightarrow$ 50% Heads, 50% Tails
2. An unbiased flip(er) means
   • each flip is independent of the others $$P(\mathcal{F}_1(\mathcal{X})*\mathcal{F}_2(\mathcal{X}))= P(\mathcal{F}_1(\mathcal{X}))*P(\mathcal{F}_2(\mathcal{X}))$$
   • the expected value of the flip is the same as that of the coin $$E(\prod_{i=0}^\infty \mathcal{F}_i(\mathcal{X})) = E(\mathcal{X})$$

Simple Model to Reality

Coin Flip

1. Each flip gives us a small piece of information about the flipper and the coin
2. More flips provides more information
   • Random / Stochastic variations in coin and flipper cancel out
   • Systematic variations accumulate

Real Experiment

1. Each measurement tells us about our sample, out instrument, and our analysis
2. More measurements provide more information
   • Random / Stochastic variations in sample, instrument, and analysis cancel out
   • Normally the analysis has very little to no stochastic variation
   • Systematic variations accumulate

Reproducibility

A very broad topic with plenty of sub-areas and deeper meanings. We mean two things by reproducibility

Analysis

The process of going from images to numbers is detailed in a clear manner that anyone, anywhere could follow and get the exact (within some tolerance) same numbers from your samples

• No platform dependence
• No proprietary or "in house" algorithms
• No manual clicking, tweaking, or copying
• One script to go from image to result

Measurement

Everything for analysis + taking a measurement several times (noise and exact alignment vary each time) does not change the statistics significantly

• No sensitivity to mounting or rotation
• No sensitivity to noise
• No dependence on exact illumination

Reproducible Analysis

The basis for reproducible scripts and analysis are scripts and macros. Since we will need to perform the same analysis many times to understand how reproducible it is.

IMAGEFILE=$1
THRESHOLD=130
matlab -r "inImage=$IMAGEFILE; threshImage=inImage>$THRESHOLD; analysisScript;"

• or java -jar ij.jar -macro TestMacro.ijm blobs.tif
• or Rscript -e "library(plyr);..."

Comparing Groups: Intraclass Correlation Coefficient

The intraclass correlation coefficient basically looking at how similar objects within a group are compared to between groups

+/- r Code

test.data<-ldply(c("Cats","People","ID","Tiger"),function(class.name) {
  data.frame(name=class.name,values=runif(40))
})
ggplot(test.data,aes(x=name,y=values))+
  geom_point()+
  labs(x="Groups",y="Value",title="Low Group Similarity")+
  theme_bw(20)

+/- r Code

test.data<-ldply(c("Cats","People","ID","Tiger"),function(class.name) {
  data.frame(name=class.name,values=nchar(class.name)+runif(40))
})
ggplot(test.data,aes(x=name,y=values))+
  geom_point()+
  labs(x="Groups",y="Value",title="High Group Similarity")+
  theme_bw(20)

Intraclass Correlation Coefficient Definition

$$ICC = \frac{S_A^2}{S_A^2+S_W^2}$$

where

• $S_A^2$ is the variance among groups or classes
  • Estimate with the standard deviations of the mean values for each group

• $S_W^2$ is the variance within groups or classes.

  • Estimate with the average of standard deviations for each group

• 1 means 100% of the variance is between classes

• 0 means 0% of the variance is between classes

Intraclass Correlation Coefficient: Values

+/- r Code

library("ICC")

test.data<-ldply(c("Cats","People","ID","Tiger"),function(class.name) {
  data.frame(name=class.name,values=runif(40))
})
icVal<-ICCbare(name,values,data=test.data)$ICC
ggplot(test.data,aes(x=name,y=values))+
  geom_point()+
  labs(x="Groups",y="Value",title=paste("Low Group Similarity\n ICC:",icVal))+
  theme_bw(20)

+/- r Code

test.data<-ldply(c("Cats","People","ID","Tiger"),function(class.name) {
  data.frame(name=class.name,values=nchar(class.name)+runif(40))
})
icVal<-ICCbare(name,values,data=test.data)$ICC
ggplot(test.data,aes(x=name,y=values))+
  geom_point()+
  labs(x="Groups",y="Value",title=paste("High Group Similarity\n ICC:",icVal))+
  theme_bw(20)

Intraclass Correlation Coefficient: Values for Coin-Flips

We have one biased coin and try to figure out how many flips we need for the ICC to tell the difference to the normal coin

+/- r Code

library("ICC")
name.list<-c("Coin A","Coin B")
test.data<-ldply(c(1:length(name.list)),function(class.id) {
  data.frame(name=name.list[class.id],values=runif(100,max=1+2*(class.id-1))>0.5)
})
icVal<-ICCbare(name,values,data=test.data)$ICC
ggplot(test.data,aes(x=name,y=values))+
  geom_jitter()+
  labs(x="Groups",y="Value",title=paste("100 flips\n ICC:",icVal))+
  theme_bw(20)

With many thousands of flips we eventually see a very strong difference but unless it is very strongly biased ICC is a poor indicator for the differences

+/- r Code

name.list<-c("Coin A","Coin B")
test.data<-ldply(c(1:length(name.list)),function(class.id) {
  data.frame(name=name.list[class.id],values=runif(20000,max=1+2*(class.id-1))>0.5)
})
icVal<-ICCbare(name,values,data=test.data)$ICC
ggplot(test.data,aes(x=name,y=values))+
  geom_jitter()+
  labs(x="Groups",y="Value",title=paste("20,000 flips\n ICC:",icVal))+
  theme_bw(20)

Comparing Groups: Tests

Once the reproducibility has been measured, it is possible to compare groups. The idea is to make a test to assess the likelihood that two groups are the same given the data

1. List assumptions
2. Establish a null hypothesis
   • Usually both groups are the same
3. Calculate the probability of the observations given the truth of the null hypothesis
   • Requires knowledge of probability distribution of the data
   • Modeling can be exceptionally complicated

Loaded Coin

We have 1 coin from a magic shop

• our assumptions are
  • we flip and observe flips of coins accurately and independently
  • the coin is invariant and always has the same expected value
• our null hypothesis is the coin is unbiased $E(\mathcal{X})=0.5$
• we can calculate the likelihood of a given observation given the number of flips (p-value)

n.flips<-c(1,5,10)
cf.table<-data.frame(No.Flips=n.flips,PAH=paste(round(1000*0.5^n.flips)/10,"%"))
names(cf.table)<-c("Number of Flips","Probability of All Heads Given Null Hypothesis (p-value)")
kable(cf.table)

How good is good enough?

Comparing Groups: Student's T Distribution

Since we do not usually know our distribution very well or have enough samples to create a sufficient probability model

Student T Distribution

We assume the distribution of our stochastic variable is normal (Gaussian) and the t-distribution provides an estimate for the mean of the underlying distribution based on few observations.

• We estimate the likelihood of our observed values assuming they are coming from random observations of a normal process

Student T-Test

Incorporates this distribution and provides an easy method for assessing the likelihood that the two given set of observations are coming from the same underlying process (null hypothesis)

• Assume unbiased observations
• Assume normal distribution

Multiple Testing Bias

Back to the magic coin, let's assume we are trying to publish a paper, we heard a p-value of < 0.05 (5%) was good enough. That means if we get 5 heads we are good!

n.flips<-c(1,4,5)
cf.table<-data.frame(No.Flips=n.flips,PAH=paste(round(1000*0.5^n.flips)/10,"%"))
names(cf.table)<-c("Number of Flips","Probability of All Heads Given Null Hypothesis (p-value)")
kable(cf.table)

n.friends<-c(1,10,20,40,80)
cfr.table<-data.frame(No.Friends=n.friends,PAH=paste(round((1000*(1-(1-0.5^5)^n.friends)))/10,"%"))
names(cfr.table)<-c("Number of Friends Flipping","Probability Someone Flips 5 heads")
kable(cfr.table)

Clearly this is not the case, otherwise we could keep flipping coins or ask all of our friends to flip until we got 5 heads and publish

The p-value is only meaningful when the experiment matches what we did.

• We didn't say the chance of getting 5 heads ever was < 5%
• We said if we have exactly 5 observations and all of them are heads the likelihood that a fair coin produced that result is <5%

Many methods to correct, most just involve scaling $p$. The likelihood of a sequence of 5 heads in a row if you perform 10 flips is 5x higher.

Multiple Testing Bias: Experiments

This is very bad news for us. We have the ability to quantify all sorts of interesting metrics

• cell distance to other cells
• cell oblateness
• cell distribution oblateness

So lets throw them all into a magical statistics algorithm and push the publish button

With our p value of less than 0.05 and a study with 10 samples in each group, how does increasing the number of variables affect our result

make.random.data<-function(n.groups=2,n.samples=10,n.vars=1,rand.fun=runif,group.off.fun=function(grp.id) 0) {
  ldply(1:n.groups,function(c.group) {
    data.frame(group=c.group,
               do.call(cbind,llply(1:n.vars,function(c.var) group.off.fun(c.group)+rand.fun(n.samples)))
               )
  })
}
# only works for two groups
all.t.test<-function(in.data) {
  group1<-subset(in.data,group==1)[,-1,drop=F]
  group2<-subset(in.data,group==2)[,-1,drop=F]
  ldply(1:ncol(group1),function(var.id) {
    tres<-t.test(group1[,var.id],group2[,var.id])
    data.frame(var.col=var.id,
              p.val=tres$p.value,
               method=tres$method,
              var.count=ncol(group1),
              sample.count=nrow(in.data))
    }
  )
}
# run the entire analysis several times to get an average
test.random.data<-function(n.test=10,...) {
  ldply(1:n.test,function(c.test) cbind(test.num=c.test,all.t.test(make.random.data(...))))
}

var.range<-round(seq(1,60,length.out=15))
test.cnt<-80
sim.data<-ldply(var.range,function(n.vars) test.random.data(test.cnt,n.vars=n.vars))
sig.likelihood<-ddply(sim.data,.(var.count),function(c.tests) {
  data.frame(sig.vars=nrow(subset(c.tests,p.val<=0.05))/length(unique(c.tests$test.num)))
})

+/- r Code

ggplot(sig.likelihood,
       aes(x=var.count,y=sig.vars))+
  geom_point()+geom_line()+
  labs(x="Number of Variables in Study",y="Number of Significant \n (P<0.05) Findings")+
  theme_bw(20)

Multiple Testing Bias: Correction

Using the simple correction factor (number of tests performed), we can make the significant findings constant again

+/- r Code

sig.likelihood.corr<-ddply(sim.data,.(var.count),function(c.tests) {
  data.frame(sig.vars=nrow(subset(c.tests,p.val<=0.05/var.count))/length(unique(c.tests$test.num)))
})
ggplot(sig.likelihood.corr,
       aes(x=var.count,y=sig.vars))+
  geom_point()+geom_line(aes(color="Corrected"))+
  geom_point(data=sig.likelihood)+
  geom_line(data=sig.likelihood,aes(color="Non-Corrected"))+
  geom_hline(yintercept=0.05,color="green",alpha=0.4,size=2)+
  scale_y_sqrt()+
  labs(x="Number of Variables in Study",y="Number of Significant \n (P<0.05) Findings")+
  theme_bw(20)

So no harm done there we just add this correction factor right? Well what if we have exactly one variable with shift of 1.0 standard deviations from the other.

var.range<-round(seq(10,60,length.out=10))
test.cnt<-100
one.diff.sample<-function(grp.id) ifelse(grp.id==2,.10,0)
sim.data.diff<-ldply(var.range,function(n.samples) 
  test.random.data(test.cnt,n.samples=n.samples,
                   rand.fun=function(n.cnt) rnorm(n.cnt,mean=1,sd=0.1),
                   group.off.fun=one.diff.sample))

+/- r Code

ggplot(sim.data.diff,aes(x=sample.count,y=p.val))+
  geom_point()+
  geom_smooth(aes(color=" 1 Variable"))+
  geom_hline(yintercept=0.05,color="green",alpha=0.4,size=2)+
  labs(x="Number of Samples in Study",y="P-Value for a 10% Difference")+
  theme_bw(20)

Multiple Testing Bias: Sample Size

var.range<-c(1,5,10,20,100) # variable count
sim.data.psig<-ldply(var.range,function(c.vcnt) {
  cbind(var.count=c.vcnt,ddply(sim.data.diff,.(sample.count),function(c.sample) 
    data.frame(prob.sig=nrow(subset(c.sample,p.val<=0.05/c.vcnt))/nrow(c.sample))
  ))
})

+/- r Code

ggplot(sim.data.psig,aes(x=sample.count,y=100*prob.sig))+
  geom_line(aes(color=as.factor(var.count)),size=2)+
  ylim(0,100)+
  labs(x="Number of Samples in Study",y="Probability of Finding\n Significant Variable (%)",color="Variables")+
  theme_bw(20)

Parameters

make.im.proc.chain<-function(root.node="Raw\nImages",filters=c(),filter.parms=c(),
                             segmentation=c(),segmentation.parms=c(),
                             analysis=c(),analysis.parms=c()) {
  node.names<-c("Raw\nImages",
                filter.parms,filters,
                segmentation.parms,segmentation,
                analysis.parms,analysis
                
                )
  
  c.mat<-matrix(0,length(node.names),length(node.names))
  colnames(c.mat)<-node.names
  rownames(c.mat)<-node.names
  
  
  for(cFilt in filters) {
    c.mat["Raw\nImages",cFilt]<-1
    for(cParm in filter.parms) c.mat[cParm,cFilt]<-1
    for(cSeg in segmentation) {
        c.mat[cFilt,cSeg]<-1
        for(cParm in segmentation.parms) c.mat[cParm,cSeg]<-1
        for(cAnal in analysis) {
          c.mat[cSeg,cAnal]<-1
          for(cParm in analysis.parms) c.mat[cParm,cAnal]<-1
        }
      }
    }
  
  
  g<-graph.adjacency(c.mat,mode="directed")
  V(g)$degree <- degree(g)
  V(g)$label <- V(g)$name
  V(g)$color <- "lightblue"
  V(g)["Raw\nImages"]$color<-"lightgreen"
  for(cAnal in analysis) V(g)[cAnal]$color<-"pink"
  V(g)$size<-30
  for(cParam in c(filter.parms,segmentation.parms,analysis.parms)) {
    V(g)[cParam]$color<-"grey"
    V(g)[cParam]$size<-25
  }
  E(g)$width<-2
  g
  }

+/- r Code

g<-make.im.proc.chain(filters=c("Gaussian\nFilter"),
                      filter.parms=c("3x3\nNeighbors","0.5 Sigma"),
                      segmentation=c("Threshold"),
                      segmentation.parms=c("100"),
                      analysis=c("Shape\nAnalysis","Thickness\nAnalysis")
                      )
plot(g)#,layout=layout.circle) #, layout=layout.circle)# layout.fruchterman.reingold)# layout.kamada.kawai) 

• Green are the images we start with (measurements)
• Blue are processing steps
• Gray are use input parameters
• Pink are the outputs

The Full Chain

+/- r Code

library(igraph)
g<-make.im.proc.chain(filters=c("Gaussian\nFilter","Median\nFilter","Diffusion\nFilter","No\nFilter",
                                "Laplacian\nFilter"),
                      segmentation=c("Threshold","Hysteresis\nThreshold","Automated"),
                      analysis=c("Shape\nAnalysis","Thickness\nAnalysis","Distribution\nAnalysis",
                                 "Skeleton\nAnalysis","2 Point\nCorr","Curvature")
                      )
plot(g,layout=layout.reingold.tilford) #, layout=layout.circle)# layout.fruchterman.reingold)# layout.kamada.kawai) 

The Full Chain (with Parameters)

+/- r Code

g<-make.im.proc.chain(filters=c("Gaussian\nFilter","Median\nFilter","Diffusion\nFilter"),
                      filter.parms=c("3x3\nNeighbors","5x5\nNeighbors","7x7\nNeighbors",
                                     "0.5 Sigma","1.0 Sigma","1.2 Sigma"),
                      segmentation=c("Threshold","Hysteresis\nThreshold","Automated"),
                      segmentation.parms=paste(seq(90,110,length.out=3)),
                      analysis=c("Shape\nAnalysis","Thickness\nAnalysis","Distribution\nAnalysis","Skeleton\nAnalysis","2 Point\nCorr")
                      )
plot(g,layout=layout.lgl(g,maxiter=10000,root=1)) #, layout=layout.circle)# layout.fruchterman.reingold)# layout.kamada.kawai) 

• A mess, over 1080 combinations for just one sample (not even exploring a very large range of threshold values)
• To calculate this for even one sample can take days (weeks, years)
  • 512 x 512 x 512 foam sample $\rightarrow$ 12 weeks of processing time
  • 1024 x 1024 x 1024 femur bone $\rightarrow$ 1.9 years
• Not all samples are the same
• Once the analysis is run we have a ton of data
  • femur bone $\rightarrow$ 60 million shapes analyzed
• What do we even want?
• How do we judge the different results?

Qualitative vs Quantitative

Given the complexity of the tree, we need to do some pruning

Qualitative Assessment

• Evaluating metrics using visual feedback
• Compare with expectations from other independent techniques or approach
• Are there artifacts which are included in the output?
• Do the shapes look correct?
• Are they distributed as expected?
• Is their orientation meaningful?

+/- r Code

title.fun<-function(file.name) ""
show.pngs.as.grid(Sys.glob("ext-figures/poros*.png"),title.fun,zoom=0.5)

Quantitative Metrics

With a quantitative approach, we can calculate the specific shape or distribution metrics on the sample with each parameter and establish the relationship between parameter and metric.

Parameter Sweep

The way we do this is usually a parameter sweep which means taking one (or more) parameters and varying them between the reasonable bounds (judged qualitatively).

source('~/Dropbox/tipl_maven/snippets/R/shapeAnalysisProcess.R')
source('~/Dropbox/tipl_maven/snippets/R/commonReportFunctions.R')
# read and correct the coordinate system
thresh.fun<-function(x) {
  pth<-rev(strsplit(x,"/")[[1]])[2]
  t<-strsplit(pth,"_")[[1]][3]
  as.numeric(substring(t,2,nchar(t)))
}
readfcn<-function(x) cbind(compare.foam.corrected(x,
                                                  checkProj=F
                                                  #force.scale=0.011 # force voxel size to be 11um
                                                  ),
                           thresh=thresh.fun(x) # how to parse the sample names
                           )
# Where are the csv files located
rootDir<-"../common/data/mcastudy" 
clpor.files<-Sys.glob(paste(rootDir,"/a*/lacun_0.csv",sep="/")) # list all of the files

# Read in all of the files
all.lacun<-ldply(clpor.files,readfcn,.parallel=T)

+/- r Code

 ggplot(all.lacun,aes(y=VOLUME*1e9,x=thresh))+
  geom_jitter(alpha=0.1)+geom_smooth()+
  theme_bw(24)+labs(y="Volume (um3)",x="Threshold Value",color="Threshold")+ylim(0,1000)

Is it always the same?

+/- r Code

 ggplot(subset(all.lacun,thresh %% 1000==0),aes(y=VOLUME*1e9,x=as.factor(thresh)))+
  geom_violin()+
  theme_bw(24)+labs(y="Volume (um3)",x="Threshold Value",color="Threshold")+ylim(0,1000)

+/- r Code

 ggplot(all.lacun,aes(y=PCA1_Z,x=thresh))+
  geom_jitter(alpha=0.1)+geom_smooth()+
  theme_bw(24)+labs(y="Orientation",x="Threshold Value",color="Threshold")

Sensitivity

Sensitivity is defined in control system theory as the change in the value of an output against the change in the input. $$S = \frac{|\Delta \textrm{Metric}|}{|\Delta \textrm{Parameter}|}$$

Such a strict definition is not particularly useful for image processing since a threshold has a unit of intensity and a metric might be volume which has $m^3$ so the sensitivity becomes volume per intensity.

Practical Sensitivity

A more common approach is to estimate the variation in this parameter between images or within a single image (automatic threshold methods can be useful for this) and define the sensitivity based on this variation. It is also common to normalize it with the mean value so the result is a percentage.

$$S = \frac{max(\textrm{Metric})-min(\textrm{Metric})}{avg(\textrm{Metric})}$$

Sensitivity: Real Measurements

In this graph it is magnitude of the slope. The steeper the slope the more the metric changes given a small change in the parameter

+/- r Code

poresum<-function(all.data) ddply(all.data,.(thresh),function(c.sample) {
  data.frame(Count=nrow(c.sample),
             Volume=mean(c.sample$VOLUME*1e9),
             Stretch=mean(c.sample$AISO),
             Oblateness=mean(c.sample$OBLATENESS),
             #Lacuna_Density_mm=1/mean(c.sample$DENSITY_CNT),
             Length=mean(c.sample$PROJ_PCA1*1000),
             Width=mean(c.sample$PROJ_PCA2*1000),
             Height=mean(c.sample$PROJ_PCA3*1000),
             Orientation=mean(abs(c.sample$PCA1_Z)))
})
comb.summary<-cbind(poresum(all.lacun),Phase="Lacuna")
splot<-ggplot(comb.summary,aes(x=thresh))
splot+geom_line(aes(y=Count))+geom_point(aes(y=Count))+scale_y_log10()+
  theme_bw(24)+labs(y="Object Count",x="Threshold",color="Phase")

Comparing Different Variables we see that the best (lowest) value for the count sensitivity is the highest for the volume and anisotropy.

+/- r Code

calc.sens<-function(in.df) {
  data.frame(sens.cnt=100*with(in.df,(max(Count)-min(Count))/mean(Count)),
             sens.vol=100*with(in.df,(max(Volume)-min(Volume))/mean(Volume)),
             sens.stretch=100*with(in.df,(max(Stretch)-min(Stretch))/mean(Stretch))
             )
}
sens.summary<-ddply.cutcols(comb.summary,.(cut_interval(thresh,5)),calc.sens)
ggplot(sens.summary,aes(x=thresh))+
  geom_line(aes(y=sens.cnt,color="Count"))+
  geom_line(aes(y=sens.vol,color="Volume"))+
  geom_line(aes(y=sens.stretch,color="Anisotropy"))+
  labs(x="Threshold",y="Sensitivity (%)",color="Metric")+
  theme_bw(20)

Which metric is more important?

Unit Testing

In computer programming, unit testing is a method by which individual units of source code, sets of one or more computer program modules together with associated control data, usage procedures, and operating procedures, are tested to determine if they are fit for use.

• Intuitively, one can view a unit as the smallest testable part of an application
• Unit testing is possible with every language
• Most (Java, C++, Matlab, R, Python) have built in support for automated testing and reporting

Unit Testing Continued

The first requirement for unit testing to work well is to have you tools divided up into small independent parts (functions)

• Each part can then be tested independently (unit testing)
  • If the tests are well done, units can be changed and tested independently
  • Makes upgrading or expanding tools easy
• The entire path can be tested (integration testing)
  • Catches mistakes in integration or glue

Ideally with realistic but simulated test data

• The utility of the testing is only as good as the tests you make

Example

• Given the following function function vxCnt=countVoxs(inImage)
• We can right the following tests
  • testEmpty2d

assert countVoxs(zeros(3,3)) == 0

• testEmpty3d

assert countVoxs(zeros(3,3,3)) == 0

• testDiag3d

assert countVoxs(eye(3)) == 3

Unit Testing: Examples

• Given the following function function shapeTable=shapeAnalysis(inImage) We should decompose the task into sub-components

• function clImage=componentLabel(inImage)

• function objInfo=analyzeObject(inObject)

  • function vxCnt=countVoxs(inObject)
  • function covMat=calculateCOV(inObject)
  • function shapeT=calcShapeT(covMat)
  • function angle=calcOrientation(shapeT)
  • function aniso=calcAnisotropy(shapeT)

Unit Testing in ImageJ

On this page

Unit Testing in KNIME

Read more and Here The Java-based unit-testing can be used (JUnit) before any of the plugins are compiled, additionally entire workflows can be made to test the objects using special testing nodes like

• difference node (check the if two values are different)

• disturber node (insert random / missing values to determine fault tolerance)

Test Driven Programming

Test Driven programming is a style or approach to programming where the tests are written before the functional code. Like very concrete specifications. It is easy to estimate how much time is left since you can automatically see how many of the tests have been passed. You and your collaborators are clear on the utility of the system.

1. shapeAnalysis must give an anisotropy of 0 when we input a sphere
2. shapeAnalysis must give the center of volume within 0.5 pixels
3. shapeAnalysis must run on a 1000x1000 image in 30 seconds

Visualization

One of the biggest problems with big sciences is trying to visualize a lot of heterogeneous data.

• Tables are difficult to interpret
• 3D Visualizations are very difficult to compare visually
• Contradictory necessity of simple single value results and all of the data to look for trends and find problems

Bad Graphs

There are too many graphs which say

• ‘my data is very complicated’
• ‘I know how to use __ toolbox in Matlab/R/Mathematica’

• Most programs by default make poor plots
  • Good visualizations takes time

Key Ideas

1. What is my message?
2. Does the graphic communicate it clearly?
3. Is a graphic representation really necessary?
4. Does every line / color serve a purpose?
   • Pretend ink is very expensive

Simple Rules

1. Never use 3D graphics when it can be avoided (unless you want to be deliberately misleading), our visual system is not well suited for comparing heights of different
2. Pie charts can also be hard to interpret
3. Background color should almost always be white (not light gray)
4. Use color palettes adapted to human visual sensitivity

What is my message

• Plots to "show the results" or "get a feeling" are usually not good

+/- r Code

xd<-runif(80)
test.data<-data.frame(x=xd,y=xd+runif(80),z=runif(80))
plot(test.data)

• Focus on a single, simple message
  • X is a little bit correlated with Y
    

    +/- r Code

    ggplot(test.data,aes(x,y))+
    geom_point()+geom_smooth(method="lm")+
    coord_equal()+
    labs(title="X is weakly correlated with Y")+
    theme_bw(20)
    

Does my graphic communicate it clearly?

• Too much data makes it very difficult to derive a clear message
  

  +/- r Code

  xd<-runif(5000)
  test.data<-data.frame(x=xd,y=(xd-0.5)*runif(5000))
  ggplot(test.data,aes(x,y))+
  geom_point()+
  coord_equal()+
  theme_bw(20)
  

• Filter and reduce information until it is extremely simple

+/- r Code

ggplot(test.data,aes(x,y))+
  stat_binhex(bins=20)+
  geom_smooth(method="lm",aes(color="Fit"))+
  coord_equal()+
  theme_bw(20)+guides(color=F)

+/- r Code

ggplot(test.data,aes(x,y))+
  geom_density2d(aes(color="Contour"))+
  geom_smooth(method="lm",aes(color="Linear Fit"))+
  coord_equal()+
  labs(color="Type")+
  theme_bw(20)

Grammar of Graphics

• What is a grammar?
  • Set of rules for constructing and validating a sentence
  • Specifies the relationship and order between the words constituting the sentence
• How does this apply to graphics?
  • If we develop a consistent way of expressing graphics (sentences) in terms of elements (words) we can compose and decompose graphics easily

The most important modern work in graphical grammars is “The Grammar of Graphics” by Wilkinson, Anand, and Grossman (2005). This work built on earlier work by Bertin (1983) and proposed a grammar that can be used to describe and construct a wide range of statistical graphics.

H. Wickham. ggplot2: elegant graphics for data analysis. Springer New York, 2009.

Grammar Explained

Normally we think of plots in terms of some sort of data which is fed into a plot command that produces a picture

• In Excel you select a range and plot-type and click "Make"
• In Matlab you run plot(xdata,ydata,color/shape)

1. These produces entire graphics (sentences) or at least phrases in one go and thus abstract away from the idea of grammar.
2. If you spoke by finding entire sentences in a book it would be very ineffective, it is much better to build up word by word

Grammar

Separate the graph into its component parts

1. Data Mapping
   • $var1 \rightarrow x$, $var2 \rightarrow y$

1. Points
2. Axes / Coordinate System
3. Labels / Annotation

Construct graphics by focusing on each portion independently.

Wrapping up

• I am not a statistician
• This is not a statistics course
• If you have questions or concerns, Both ETHZ and Uni Zurich offer free consultation with real statisticians
  • They are rarely bearers of good news

• Simulations are very helpful
• Try and understand the tests you are performing