## Power Analysis and the Probability of Errors

Power analysis is a very useful tool to estimate the statistical power from a study. It effectively allows a researcher to determine the needed sample size in order to obtained the required statistical power. Clients often ask (and rightfully so) what the sample size should be for a proposed project. Sample sizes end up being a delicate balance between the amount of acceptable error, detectable effect size, power, and financial cost. A lot of factors go into this decision. This example will discuss one approach.

In more specific terms power is the probability that a statistical test will reject the null hypothesis when the null hypothesis is truly false. What this means is that when power increases the probability of making a Type II error decreases. The probability of a Type II error is denoted by $\beta$ and power is calculated as $power = 1 - \beta$.

In order to calculate the probability of a Type II error a researcher needs to know a few pieces of information $\mu$, $\sigma^2$, $n$, and $\alpha$ (probability of a Type I error). Normally, if a researcher already knows the population mean ($\mu$) and variance ($\sigma^2$) there is no need to take a sample to estimate them. However, we can set it up so we can look at a range of possible unknown population means and variances to see what the probability of a Type II error is for those values.

The following code shows a basic calculation and the density plot of a Type II error.

$P\left(X>38| \mu,\sigma^2 \right)=P\left(Z>\frac{38-50}{10}\right)=P\left(Z>-1.2\right)=.8849303$

my.means = c(38,50);
my.sd = c(10,10);

X = seq(my.means[1]-my.sd[1]*5,my.means[1]+my.sd[1]*5, length=100);
dX = dnorm(X, mean=my.means[1], sd=my.sd[1]); #True distribution
dX = dX/max(dX);

X.2 = seq(my.means[2]-my.sd[2]*5,my.means[2]+my.sd[2]*5, length=100);
dX.2 = dnorm(X.2, mean=my.means[2], sd=my.sd[2]); #Sampled distribution
dX.2 = dX.2/max(dX.2);

plot(X, dX, type="l", lty=1, lwd=2, xlab="x value",
ylab="Density", main="Comparison of Type II error probability",
xlim=c(min(X,X.2),max(X,X.2)))
lines(X.2,dX.2, lty=2)
abline(v=my.means, lty=c(1,2));
legend("topleft", c("Sample","Type II Error"), lty=c(2,1), col=c("black","black"), lwd=2, title=expression(paste("Type II Error")))

x = (my.means[1]-my.means[2])/my.sd[2];
1-pnorm(x);

[1] 0.8849303



However, there is a more direct way to determine necessary power and the Type II error calculation using R by setting up a range of sample sizes. Often the mean and variance for power analysis are established through a pilot study or other preliminary research on the population of interest.

R has several convenient functions to calculate power for many different types of tests and designs.  For this example power.anova.test will be used for demonstration purposes.  However, power.t.test could have been used as well.

For each of these functions there are varying values that need to be supplied however they generally involve effect size (difference in group means), sample size, alpha level, measure of variability and power.  In this example the within variability is estimated by using the MSE from the ANOVA table and the between variability is estimated from the variance of group means.  Otherwise, three different alpha levels are measured over a range of 50 possible n’s.


set.seed(1234);
x = matrix(NA,nrow=200,ncol=2);
x[,1]= rbind(rep(1,100),rep(2,100))
x[,2]= rbind(rnorm(100,23,20),rnorm(100,7,20));
x=data.frame(x); names(x) = c("group","value");
groupmeans = as.matrix(by(x$value,x$group,mean));
x.aov = summary(aov(x$value ~ as.factor(x$group)))

nlen = 50
withinvar = 470; ##From MSE in ANOVA
raw.pwr = matrix(NA, nrow=nlen, ncol=5)
for(i in 1:nlen){
pwr.i10 = power.anova.test(groups=length(groupmeans), n=1+i,
between.var=var(groupmeans), within.var=withinvar, sig.level=.10)
pwr.i05 = power.anova.test(groups=length(groupmeans), n=1+i,
between.var=var(groupmeans), within.var=withinvar, sig.level=.05)
pwr.i01 = power.anova.test(groups=length(groupmeans), n=1+i,
between.var=var(groupmeans), within.var=withinvar, sig.level=.01)
power10 = pwr.i10$power power05 = pwr.i05$power

## Hurricane Sandy Land Wind Speed and Kriging

NJ Hurricane Sandy Landfall Data

These data come from the National Climatic Data Center (NCDC).  Using the above link will download all of the data collected by the NCDC on the day of Hurricane Sandy.  The data can also be obtained directly from the source at http://cdo.ncdc.noaa.gov/qclcd/QCLCD.

The purpose of this post is not a discussion on kriging or any of its properties. The purpose is to simply provide a simple R example on kriging and how it can be applied on real data.   This R code uses Hurricane Sandy as an example and predicts the wind speed across the state of New Jersey using the measurement stations that the NCDC uses. It is easy to see how this technique can be applied to other data.

wd <- "C:\\nj hurricane data" setwd(wd) my.files <- dir(wd, pattern = ".txt", full.names = TRUE, ignore.case = TRUE) krig <- matrix(NA, nrow=length(my.files), ncol=3) for(i in 1:length(my.files)){ raw <- read.csv(my.files[i], skip=6, header=T) raw <- subset(raw, raw$Date=="20121029") latlongraw <- read.csv(my.files[i], skip=3, header=F, sep="") latlongrows <- head(latlongraw,2)[,2] latlong <- as.numeric(levels(latlongrows)[latlongrows]) raw$WindSpeed[raw$WindSpeed=="M"] <- NA if(is.factor(raw$WindSpeed)){ raw$WindSpeed <- as.numeric(levels(raw$WindSpeed)[raw$WindSpeed]) } krig[i,1] <- latlong[1] krig[i,2] <- latlong[2] krig[i,3] <- max(raw$WindSpeed, na.rm=T) } lat <- krig[,1] long <- krig[,2] s<-cbind(long,lat) PM <- krig[,3] sandy$coords <- s sandy$data <- PM ml <- likfit(sandy, fix.nugget=F, cov.model="exponential", ini = c(10, 5), nugget=4) ml summary(ml) grid <- map("state","new jersey", plot=FALSE) grid.plot <- cbind( c(min(grid$x, na.rm=TRUE),max(grid$x, na.rm=TRUE)), c(min(grid$y, na.rm=TRUE),max(grid$y, na.rm=TRUE)) ) sp1<-seq(grid.plot[1,1]-.25,grid.plot[2,1]+.25,length=100) sp2<-seq(grid.plot[1,2]-.25,grid.plot[2,2]+.25,length=100) sp<-expand.grid(sp1,sp2) inLoc<-map.where("state",x=sp[,1],y=sp[,2]) inLoc[is.na(inLoc)]<-"NA" inLoc<-inLoc=="new jersey" #Perform ordinary Kriging (value of cov.pars and nugget are copied from mle output): pred<-krige.conv(data=PM,coords=s,locations=sp, krige=krige.control(type.krige="ok",cov.model="exponential", cov.pars=c(100,3), nugget=5)) pred$predict[!inLoc]<-NA pred$krige.var[!inLoc]<-NA #Plot the predicted values: image.plot(sp1,sp2,matrix(pred$predict,100,100),zlim=range(PM), main="Sandy Maximum Wind Speed in MPH", xlab="Longitude",ylab="Latitude") map("county",add=T) points(s, pch=16) [/sourcecode] ## What Time Is It? A common scenario that I run into is time and how to deal with it. I often will do a variety of summaries and analysis that need to be measured at different points in time. Whether I want to graph the data or review the results I need to be able to perform measurements relative to time and interpret the time output in human readable form. R has several functions to handle those time related scenario. For starters one thing needs to be acknowledged. On Unix based system the beginning of time is January 1, 1970. At that moment exactly zero (0) seconds have passed and it is known as Unix epoch. Anything before that is negative and anything after that time is positive. This measure of time is defined in seconds. This piece of code shows the Unix epoch and how it compares to the current time. In this example the origin is zero (0) but as one can see it can easily be changed to calculate the length of time that has past since any given time. ##Local Time in EDT (GMT-5). If local time IS GMT then the date would be January 1, 1970 iso <- ISOdatetime(1969,12,31,19,0,0) # (YYYY,MM,DD,HH,MM,SS) origin <- as.double(iso) t <- Sys.time() curr.time <- as.double(t) curr.time - origin [/sourcecode] The following bit of code simulates a cumulative process where each consecutive process follows a Poisson distribution and $\lambda$ (lambda) is a random variable distributed uniformly. This would be comparable to a process where election vote counts are being collected throughout Election Night. [sourcecode language="css"] set.seed(1234) time.series <- seq(Sys.time(), Sys.time()+4*60*60, by=60*5) #Add 4 hours to time series vals <- matrix(NA, ncol=3, nrow=length(time.series)) for(i in 1:length(time.series)){ if( is.numeric(vals[i-1,1]) & is.numeric(vals[i-1,2]) ){ vals[i,1] <- vals[i-1,1]+rpois(n=1, runif(1, 700,1000)) vals[i,2] <- vals[i-1,2]+rpois(n=1, runif(1, 850,1000)) } else { vals[i,1] <- runif(1, 700,1000) vals[i,2] <- runif(1, 850,1000) } } vals[,3] <- vals[,1]/(vals[,1]+vals[,2]) prop <- vals[,3]*100 prop2 <- (1-vals[,3])*100 plot(time.series,prop, type="o", cex=.6, lty=1, col="blue", pch=1, ylim=c(46,54) , ylab="Percent", xlab="Time of Computation" , main=paste("Example of Cumulative Process Using Time") , sub="Cumulative Proportional Process", cex.sub=.60 , xaxt = "n", yaxt = "n"); lines(time.series, prop2, type="o", cex=.6, lty=1, col="red", pch=4); xaxis.seq.half<- seq(strptime(c(min(time.series)), "%Y-%m-%d %H"), strptime(c(max(time.series)), "%Y-%m-%d %H")+3600, by=1800) xaxis.seq.hour <- seq(strptime(c(min(time.series)), "%Y-%m-%d %H"), strptime(c(max(time.series)), "%Y-%m-%d %H")+3600, by=3600) axis(1, at=(xaxis.seq.hour), tcl = -0.7, lty = 1, lwd = 0.8, labels=FALSE) axis(1, at=(xaxis.seq.half), tcl = -0.3, lty = 2, lwd = 0.5, labels=FALSE) axis(2, at=seq(46,54,by=1), labels=FALSE, tcl = -0.2, lty = 1, lwd = 0.5) axis(2, at=seq(46,54,by=1), labels=TRUE, tcl = -0.7, lty = 1) axis.POSIXct(1, as.POSIXlt(xaxis.seq.hour), at=as.POSIXlt(xaxis.seq.hour), format="%H:%M", tcl = 0.3, las=0, lty = 2, lwd = 0.5, cex.axis=.7, labels=TRUE) [/sourcecode] ## One Big Number That’s Not Big Enough 32-bit systems pose an interesting problem because of some of its limitations. These systems are only capable of handling an integer that is $2^{63}-1 = 2147483647$. That is a fairly large number but it just not large enough. That is particularly true if that number is a phone number (e.g. (214)748-3647). As you can see there would be a large block of numbers that would cause problems. In other words a phone number is not an integer and should not be cast as an integer, ever! Normally on 32-bit systems if an integer exceeds that number it will either kindly return an error or it will send you back in time to the year 1901. Neither one of those options are very good. The solution: upgrade to a 64-bit system. That will allow for a maximum integer of 9,223,372,036,854,775,807. A 64-bit system will be good for another 292 million years. The time solution is then solved. However, it’s not out of the question that there could be a scenario when a researcher will exceed 9.2 quintillion. This could be the case when dealing with large data mining and trying to analyze something like each grain of sand on a beach. But I guess then it’s time to upgrade again to something bigger. Note that when dealing with numbers and computers they normally begin at 0. So at index 31 you have actually reach the 32nd position and therefore a 32-bit system. Likewise, $2^{63}-1 = 9223372036854775807$. If you’re not certain then you can quickly find out what the maximum allowable integer by using any of the following lines of code. In R you can identify the maximum integer: .Machine$integer.max


In Perl you can identify dates that go back on time:

#!/usr/local/bin/perl

use POSIX;
$ENV{‘TZ’} = “GMT”; for ($clock = 2147483641; $clock < 2147483651;$clock++) { print ctime($clock); } [/sourcecode] In PHP this is the same as the Perl script: [sourcecode language="css"] ## Using R to Compare Hurricane Sandy and Hurricane Irene Having just lived through two back to back hurricanes (Irene in 2011 and Sandy in 2012) that passed through the New York metro area I was curious how the paths of the hurricanes differed. I worked up a quick graph in R using data from Unisys. The data also includes wind speed and barometric pressure. library(maps) library(maptools) library(rgdal) sandy = read.table(file=”http://weather.unisys.com/hurricane/atlantic/2012/SANDY/track.dat”, skip=3,fill=TRUE) irene = read.table(file=”http://weather.unisys.com/hurricane/atlantic/2011H/IRENE/track.dat”,skip=3,fill=TRUE) colnames(sandy) = c(“Advisory”,”Latitude”,”Longitude”,”Time”,”WindSpeed”,”Pressure”,”Status”) colnames(irene) = c(“Advisory”,”Latitude”,”Longitude”,”Time”,”WindSpeed”,”Pressure”,”Status”) sandy$WindSpeedColor <- 'blue' sandy$WindSpeedColor[sandy$WindSpeed >= 75] <- 'red' irene$WindSpeedColor <- 'blue' irene$WindSpeedColor[sandy$WindSpeed >= 75] <- 'red' xlim <- c(-88,-65) ylim <- c(25,48) state.list <- c('new york','new jersey','virginia','massachusetts','connecticut','delaware','pennsylvania','maryland','north carolina','south carolina','georgia','florida', 'new hampshire','maine','district of columbia','west virginia','vermont') my.map <- map("state", region = state.list, interior = FALSE, xlim=xlim, ylim=ylim) map("state", region = state.list, boundary = FALSE, col="gray", add = TRUE,xlim=xlim) lines(x=sandy$Longitude,y=sandy$Latitude,col="black",cex=0.75) points(x=sandy$Longitude,y=sandy$Latitude,col=sandy$WindSpeedColor,pch=15,cex=0.9) text(x=sandy$Longitude,y=sandy$Latitude,col='dark green',labels=sandy$Pressure,adj=c(-0.9),cex=0.5) lines(x=irene$Longitude,y=irene$Latitude,col="black",cex=0.75) points(x=irene$Longitude,y=irene$Latitude,col=irene$WindSpeedColor,pch=15,cex=0.9) text(x=irene$Longitude,y=irene$Latitude,col='light green',labels=irene$Pressure,adj=c(-0.9),cex=0.5) title("Path of Hurricane Sandy (2012) and Hurricane Irene (2011)\nwith Wind Speed and Barometric Pressure") legend('topleft',c('Tropical Storm Wind Speeds','Hurricane Wind Speeds'),pch=15, col=c('blue','red')) box() [/sourcecode] ## Mapping Capabilities in R From time-to-time creating a basic map of the United States or other parts of the world to complement some statistical analysis is useful to emphasize a point. The maps package in R provide a good way to produce these these maps. These maps axes are based on latitude and longitude so overlaying other information on these maps is quite simple. Furthermore, it makes a nice addition to geo-spatial analysis. This is a simple introductory example that is designed to show a couple of different examples. First, it shows the USA lower 48 states database and how it can be used to create the state boundaries. Second, the world database is used to add Alaska and Hawaii. A nice feature is that the colors for each of the state boundaries can be easily changed to group states. This (non-optimized) code shows a basic example on how to use the maps package to show the Election Day poll closing times for all 50 states and the District of Columbia. library(maps) m <- map("state", interior = FALSE, plot=FALSE) m$my.colors <- 0 m$my.colors[m$names %in% c("virginia:main","georgia","vermont","kentucky","indiana","south carolina")] <- 2 m$my.colors[m$names %in% c("west virginia","north carolina:main","ohio")] <- 3 m$my.colors[m$names %in% c("tennessee","rhode island","oklahoma","new jersey","mississippi","massachusetts:main","illinois","maryland","maine","missouri","new hampshire","new jersey","pennsylvania","florida","connecticut","delaware","district of columbia","alabama")] <- 4 m$my.colors[m$names %in% c("arkansas")] <- 5 m$my.colors[m$names %in% c("arizona","colorado","kansas","louisiana","minnesota","michigan:north","michigan:south","nebraska","new york:main","new york:manhattan","new york:staten island","new york:long island","new mexico","north dakota","south dakota","texas","wyoming","wisconsin")] <- 6 m$my.colors[m$names %in% c("iowa","montana","nevada","utah")] <- 7 m$my.colors[m$names %in% c("california","idaho","oregon","washington:main")] <- 8 m.world <- map("world", c("USA","hawaii"), xlim=c(-180,-65), ylim=c(19,72),interior = FALSE) title("Election Day Poll Closing Times") map("state", boundary = FALSE, col="grey", add = TRUE, fill=FALSE) map("state", boundary = TRUE, col=m$my.colors, add = TRUE, fill=TRUE ) map("world", c("hawaii"), boundary = TRUE, col=8, add = TRUE, fill=TRUE ) map("world", c("USA:Alaska"), boundary = TRUE, col='orange', add = TRUE, fill=TRUE ) legend("topright", c('7:00pm','7:30pm','8:00pm','8:30pm','9:00pm','10:00pm','11:00pm','1:00am'), pch=15, col=c(2,3,4,5,6,7,8,'orange'), title="Poll Closing Time (EST)", ncol=2, cex=1.2) [/sourcecode] Because the maps are based on longitude and latitude this package can be extended beyond just county, state and country boundaries and can be used on techniques like kriging and co-kriging to develop model-based maps. ## Text Mining When it comes down to it R does a really good job handling structured data like matrices and data frames. However, its ability to work with unstructured data is still a work in progress. It can and it does handle text mining but the documentation is incomplete and the capabilities still don’t compare to other programs such as MALLET or Mahout. Though the formal documentation is still lacking. Though this is not an example on real data it does provide the basic tools on text mining and, in particular, latent dirichlet allocation. There are three R libraries that are useful for text mining: tm, RTextTools, and topicmodels. The tm library is the core of text mining capabilities in R. Unstructured text files can come in many different formats. I often find that I must get my own data and consequently the data generally originates as plain text (.txt) files. However, those who want to analyze Twitter feeds can user the twitteR library which is useful for analyzing social media topics in real time. This example will incorporate the CNN twitter feed. In order for R to interpret and analyze these text files they must ultimately be converted into a document term matrix. But first a corpus must be created. A corpus is simply a collection of documents where each document its a topic. When reading text documents directly from local file the following R code can be used. Data Preparation using Local Text Files  #These files can be just raw text. For example it could be simply copied and pasted from a Web site. dir = "C:\\Documents and Settings\\clints\\My Documents\\LDA-S"; filenames = list.files(path=dir,pattern="\\.txt"); setwd(dir); docs = NULL; titles = NULL; for (filename in filenames){ #here I specify a file that contains all the titles of the documents if(filename=="titles.txt"){ titles = paste(readLines(file(filename))); } else { docs = c(docs,list( paste(readLines(file(filename)), collapse="\n") )); } }  To pull the text from a Twitter Feed rather than text files then the following lines of code can be used. Data Preparation using Twitter  library(tm); library(RTextTools); library(topicmodels); library(twitteR); twitter_feed <- searchTwitter('@cnn', n=150); ### Optional twitter feed retrieval ##twitter_feed <- userTimeline("rdatamining", n=150); ### df <- do.call("rbind", lapply(twitter_feed, as.data.frame)); myCorpus <- Corpus(VectorSource(df$text));

k = length(docs3);
myCorpus = Corpus(VectorSource(docs));
myCorpus = tm_map(myCorpus, tolower);
myCorpus = tm_map(myCorpus, removePunctuation);
myCorpus = tm_map(myCorpus, removeNumbers);
myStopwords = c(stopwords('english'), "available", "via");
idx = which(myStopwords == "r");
myStopwords = myStopwords&#91;-idx&#93;;
myCorpus = tm_map(myCorpus, removeWords, myStopwords);

dictCorpus = myCorpus;

myCorpus = tm_map(myCorpus, stemDocument);

myCorpus = tm_map(myCorpus, stemCompletion, dictionary=dictCorpus);

myDtm = DocumentTermMatrix(myCorpus, control = list(minWordLength = 3));

findFreqTerms(myDtm, lowfreq=50);
#find the probability a word is associated
findAssocs(myDtm, 'find_a_word', 0.5);

&#91;/sourcecode&#93;

<strong>Word Cloud</strong>

library(wordcloud);
m = as.matrix(myDtm);
v = sort(colSums(m), decreasing=TRUE);
myNames = names(v);
k = which(names(v)=="miners");
myNames[k] = "mining";
d = data.frame(word=myNames, freq=v);
wordcloud(d$word, colors=c(3,4), random.color=FALSE, d$freq, min.freq=20);



Latent Dirichlet Allocation


k = 2;
SEED = 1234;
my_TM =
list(VEM = LDA(myDtm, k = k, control = list(seed = SEED)),
VEM_fixed = LDA(myDtm, k = k,
control = list(estimate.alpha = FALSE, seed = SEED)),
Gibbs = LDA(myDtm, k = k, method = "Gibbs",
control = list(seed = SEED, burnin = 1000,
thin = 100, iter = 1000)),
CTM = CTM(myDtm, k = k,
control = list(seed = SEED,
var = list(tol = 10^-4), em = list(tol = 10^-3))));

Topic = topics(my_TM[["VEM"]], 1);

#top 5 terms for each topic in LDA
Terms = terms(my_TM[["VEM"]], 5);
Terms;

(my_topics =
topics(my_TM[["VEM"]]));

most_frequent = which.max(tabulate(my_topics));

terms(my_TM[["VEM"]], 10)[, most_frequent];



Here, a model is fit setting the number of unobserved latent topics equal to two (k=2). We can then identify the most frequently occurring topics and then identify the top five terms used for the topic. In this example these are the top five terms when setting the number of groups equal to two.

 Topic 1 Topic 2 “amp” “cnn” “cnn” “tweet” “jobs” “abc” “romney” “bainport” “sensata” “cbs”

## Topic Modeling the October 2012 LDS General Conference

Twice a year, once in May, and once in October, there is the obligatory discussion on the theme of the semi-annual LDS general conference.  The importance of a conference of this type is that each person can take away key principles that will help improve their lives and that is the theme for them.

But I wanted to look at it in a more objective way for the theme of the conference.  What was the overarching theme and topics of the general conference.  There are many topics discussed but there is one that  appears more often.  After a little bit of text mining the most likely topic-terms of the October 2012 LDS General Conference is:

Christ, God, Jesus, Lord, Father, Lives

It  is not the overwhelming theme of the conference but it is often referenced.  Other topics that frequently appear include missionaries and children.  I took every word from all talks during the conference (excluding the statistical report, Priesthood, and Relief Society sessions) and analyzed the data.  I took two approaches.  The first was a simple straight forward summary of word frequency.  This ends up becoming a word cloud.  Though, personally, it’s not my favorite graph to present data the word cloud  tends to be fan-favorite and, if anything, it is a fun way to visualize data.

The second approach is solidly statistical based and uses Latent Dirichlet Allocation.  This is effectively a way to cluster the text spoken during the conference into a bunch of bag of words (that is a technical term).  Those bags contain specified words with a probability.  For those familiar with the Dirichlet distribution you will note that it is the conjugate prior of the categorical distribution.  This means that I can group the words spoken into topics (themes) and each of the groups will contain specified terms (words).

This approach shows that there is a wide distribution of topics that spanned the conference indicating that there wasn’t one dominating theme but a lot of separate themes.  However, if I were to group conference into two themes it would have to be

1) Faith in Jesus Christ

2) Time with Children.

I arrived at this conclusion by clustering the text into two groups and then based on on Latent Dirichlet Allocation found that based on that model terms relating to Christ were in one group and terms relating to children were in the other group.   Specifically the following table shows the top five terms for each topic

 Topic 1 Topic 2 Christ Children Faith Time Jesus Help God Lives Love Love

So there it is.  This is my quick analysis of the October 2012 General Conference and the theme of the conference.  However, I think the one thing that people will remember from this conference is missionary work.

*Note:  For those that care all analyses used R.