Kernal Density Plots

# When viewing this post the code looks messy, but just copy and paste into
# notepad or the R scripting window and it should line up

# This program is a little more complicated than the previous ones,and some of
# the topics may be a little out of the scope of this class, but if you
# are interested and don't understand something let me know


# example plots from Part 2 of the program using sample salary data:

# histogram of salary data with gaussian KDE plot




# KDE plots of salary data with varying bandwidth selections using the 'adjust'
# option







# PROGRAM NAME: KDE_R
# DATE: 6/18/2010
# AUTHOR : MATT BOGARD
# PURPOSE: EXAMPLE OF KERNAL DENSITY ESTIMATES IN R
# DATA USED: SIMULATED WITHIN SCRIPT
# PACKAGES: lattice Hmisc
# OUTPUT FILE: "C:\Documents and Settings\wkuuser\My Documents\TOOLS AND REFERENCES\R References"
# UPDATES:
#
#
# COMMENTS: GIVE A BASIC UNDERSTANDING OF THE density() FUNCTION IN R
# AND BANDWIDTH, KERNAL SELECTION, AND PLOTTING
#
# REFERENCES: ADAPTED FROM
# http://ise.tamu.edu/metrology/Group%20Seminars/Kernel_Density_Estimation_2008JUN27_2.pdf
#
# CONTENTS: PART1 EXAMPLES USING CODE GENERATED NORMALLY DISTRIBUTED DATA
# PART2 EXAMPLES USING SAMPLE DATA

#---------------------------------------------------------------#
# PART1 EXAMPLES USING CODE GENERATED NORMALLY DISTRIBUTED DATA
#---------------------------------------------------------------#

x <- seq(-4, 4, 0.01) y <- dnorm(x, mean=0, sd=1 # create x,y data set that is exactly normal

rn <- rnorm(1000, mean=0, sd=1 # create a set of random normal values

# plot the density of the random numbers you just generated
# let the bandwidth = nrd0 "rule of thub" h= .9*A*n*exp(-1.5), A = min{s,IQR/1.34}
# select the gaussian kernal using 1000 points (which you just generated)



plot(density(x=rn, bw="nrd0", kernel="gaussian", n=1000), col=5) # gaussian is default

lines(x, y) # plot the exactly normal data set you created earlier to compare



# plot a series of KDE plots with varying bandwitdth using 'width='option
# When this option is specified instead of 'bw',and is a number, the density
# function bandwidth is set to a kernel dependent multiple of 'width'

for(w in 1:3)
lines(density(x=rn, width=w, kernel="gaussian", n=1000), col=w+1)


# what if we change the kernal
# first close the graphics window to start over
# look at varius bandwidths using a kernal other than gaussian

# KDE otpions:
# gaussian rectangular triangular
# biweight cosine optcosine"
# epanechnikov


plot(density(x=rn, bw="nrd0", kernel="rectangular", n=1000), col=1)
for(w in 1:3)
lines(density(x=rn, width=w, kernel="rectangular", n=1000), col=w+1)


# or to just look at various densities in one plot
# using a lines statment specifying the specific kernal

plot(density(x=rn, bw="nrd0", kernel="gaussian", n=1000), col=1)
lines(density(x=rn, bw="nrd0", kernel="rectangular", n=1000), col=2)
lines(density(x=rn, bw="nrd0", kernel="triangular", n=1000), col=3)
lines(density(x=rn, bw="nrd0", kernel="epanechnikov", n=1000), col=4)
lines(density(x=rn, bw="nrd0", kernel="biweight", n=1000), col=5)
lines(density(x=rn, bw="nrd0", kernel="cosine", n=1000), col=6)
lines(density(x=rn, bw="nrd0", kernel="optcosine", n=1000), col=7)

# or to just look at various densities in one plot
# using a lines statment specifying the specific kernal

plot(density(x=rn, bw="nrd0", kernel="gaussian", n=1000), col=1)

# note if you want the default gaussian, it does not have to be specified

plot(density(x=rn, bw="nrd0", n=1000), col=1)


#-----------------------------------#
# PART2 EXAMPLES USING SAMPLE DATA
#-----------------------------------#

# load salary data

salary<-c(240,240,240,240,240,240,240,240,255,255,265,265,280,280,290,300,305,325,330,340)
print(salary)

library(Hmisc) # load for describe function
describe(salary) # general stats

library(lattice) # load lattice package for histogram graphics
hs <- hist(salary) # plot histogram and store data
print(hs)

# create a kde of salary data (gaussian)and store it as dens

dens <- density(x=salary, bw="nrd0", n=20)
print(dens)


rs <- max(hs$counts/max(dens$y)) # create scaling factor for plotting the density
lines(dens$x, dens$y*rs, col=2) # plot density over histogram using lines statment

# compare this to just plotting the non-rescaled density over the histogram

lines(density(x=salary, bw="nrd0", n=20),col=3) # green line barely shows up

# also, compare to just plotting the estimated density for the sample
# data in a separate plot


plot(density(x=salary, bw="nrd0", n=20)) # or equivalently
plot(dens, col=3)

# now that we understand the various ways we can plot a specified KDE
# we can apply the tools used in part 1 to look at the effects of
# various bandwidths and kernal choices - first close the graphing window

# what if we change the kernal
# first close the graphics window to start over
# look at varius KDEs

# KDE otpions:
# gaussian rectangular triangular
# biweight cosine optcosine"
# epanechnikov


# look at various densities in one plot
# using a lines statment specifying the specific kernal

plot(density(x=salary, bw="nrd0", kernel="gaussian", n=20), col=1)
lines(density(x=salary, bw="nrd0", kernel="rectangular", n=20), col=2)
lines(density(x=salary, bw="nrd0", kernel="triangular", n=20), col=3)
lines(density(x=salary, bw="nrd0", kernel="epanechnikov", n=20), col=4)
lines(density(x=salary, bw="nrd0", kernel="biweight", n=20), col=5)
lines(density(x=salary, bw="nrd0", kernel="cosine", n=20), col=6)
lines(density(x=salary, bw="nrd0", kernel="optcosine", n=20), col=7)

# what if we don't specify a bandwidth (clear graphics window)
# what is the default density

plot(density(x=salary, bw="nrd0", kernel="gaussian", n=20), col=1)
lines(density(x=salary, kernel="gaussian", n=20), col=2)# nrd0 appears to be default

# what if we change the default
# nrd (a variation of nrd0)
# ucv (unbiased cross validation)
# bcv (cross validation)

lines(density(x=salary, bw="ucv", kernel="gaussian", n=20), col=2)
lines(density(x=salary, bw="bcv", kernel="gaussian", n=20), col=3)
lines(density(x=salary, bw="nrd", kernel="gaussian", n=20), col=4)

# use print option to see how the bw value changes with each method

print(density(x=salary, bw="nrd0", kernel="gaussian", n=20))

# look at the differences in calculated
# bandwidth and plots with the adjust option

for(w in 1:3)
print(density(x=salary, bw="nrd0", adjust=w, kernel="gaussian", n=20))


plot(density(x=salary, bw="nrd0",adjust=1, kernel="gaussian", n=20), col=1)
lines(density(x=salary, bw="nrd0",adjust=2,kernel="gaussian", n=20), col=2)
lines(density(x=salary, bw="nrd0",adjust=3,kernel="gaussian", n=20), col=3)

Inference, Hypothesis Testing, & Regression

Inference: the exercise of providing information about the population based on information contained in the sample

Type 1 Error: rejecting a true null hypothesis
- the probabability of rejecting a true null hypothesis ( or the probability of a type 1 error) is equal to alpha or 'α' This is also the same as the level of significance in a hypothesis test.

Type II Error: failing to reject, or loosely speaking, 'accepting' the null hypthesis when it is false
-the probability of a type II error =beta or β

Power: the probability of rejecting the null hypothesis when it is false. (1-β)

Regression Model: describes how the dependent variable (y) is related to the independent variable (x), also known as the least squares line.

y = b0 + b1*x

This is derived by the process of least squares, which gives the value for the slope (b1) and the intercept (bo) that minimizes the sum of the squared deviations between the observed values of the dependent variable (Y) and the estimated values of the dependent variable (Y*).

Co-efficient of determination(R-squared) = SSR/SST -> gives the proportion of total variation explained by the regression model. Larger values indicate that the sample data are closer to the least squares line, or a stronger linear relationship exists.

Standard Normal Distribution & Z-values

############################################
# Below is R code for plotting normal distributions - illustrating
# the changes in the distribution when standard deviation changes.
# -selected graphics are imbedded below
#
# Also calculates area under the normal curve for given Z
#
# Finally plots 3-D bivariate normal density
#
#############################################



# PROGRAM NAME: NORMAL_PLOT_R
#
# ECON 206
#
# ORIGINAL SOURCE: The Standard Normal Distribution in R:
# http://msenux.redwoods.edu/math/R/StandardNormal.php
#
#
#
#


#PLOT NORMAL DENSISTY

x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=1)
plot(x,y,type="l",lwd=2,col="red")

#######################################





###################################





# INCREASE THE STANDARD DEVIATION

x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=2.5)
plot(x,y,type="l",lwd=2,col="red")


#####################################



####################################

# DECREASE THE STANDARD DEVIATION

x=seq(-4,4,length=200)
y=dnorm(x,mean=0,sd=.5)
plot(x,y,type="l",lwd=2,col="red")

# CALCULATING AREA FOR GIVEN Z-VALUES

# RECALL, FOR THE STANDARD NORMAL DISTRIBUTION THE MEAN = 0 AND
# THE STANDARD DEVIATION = 1
# THE pnorm FUNCTION GIVES THE PROBABILITY FOR THE AREA TO THE LEFT
# OF THE SPECIFIED Z-VALUE ( THE FIRST VALUE ENTERED IN THE FUNCTION)
# THE OUTPUT SHOULD MATCH WHAT YOU GET FROM THE NORMAL TABLE IN YOUR BOOK
# OR THE HANDOUT I SENT YOU

pnorm(0,mean=0, sd=1) # Z =0

pnorm(1,mean=0, sd=1) # Z <= 1

pnorm(1.55, mean=0, sd=1) # Z<=1.55

pnorm(1.645, mean=0, sd=1) # Z<= 1.645


# LETS LOOK AT A BIVARIATE NORMAL
# DISTRIBUTION


# first simulate a bivariate normal sample

library(MASS)

bivn <- mvrnorm(1000, mu = c(0, 0), Sigma = matrix(c(1, 0, 0, 1), 2))

# now we do a kernel density estimate

bivn.kde <- kde2d(bivn[,1], bivn[,2], n = 100)

# now plot your results

contour(bivn.kde)
image(bivn.kde)

persp(bivn.kde, phi = 45, theta = 30)


# fancy contour with image

image(bivn.kde); contour(bivn.kde, add = T)

# fancy perspective

persp(bivn.kde, phi = 45, theta = 30, shade = .1, border = NA)

###############################################

Regression Demo

Basic Regression in R

#
# COMMENTS: BASIC INTRODUCTION TO REGRESSION USING R
#


#
# WILLIAMS P 572 #1
#


# GET DATA

x<-c(1,2,3,4,5)

y<-c(3,7,5,11,14)


# PLOT DATA

plot(x,y)

reg1 <- lm(y~x) # COMPUTE REGRESSION ESTIMATES

summary(reg1) # PRINT OUTPUT

abline(reg1) # PLOT REGRESSION LINE



#
# SCHAUM'S P. 274 #14.2
#


#GET DATA

x<-c(20,16,34,23,27,32,18,22)
y<-c(64,61,84,70,88,92,72,77)

plot(x,y) #PLOT DATA

reg1 <- lm(y~x) # COMPUTE REGRESSION ESTIMATES

summary(reg1) # PRINT OUTPUT

abline(reg1) # PLOT REGRESSION LINE

Social Network Analysis of Tweets Using R

This past week I attended the AACRAO conference in New Orleans. This provided a great opportunity to demonstrate one of the many things you can do wit R. ( for more about the R statistical programming language see here) R is capable of many outside the box applications (beyond just the basic statistical techniques) such as social network analysis.

The image below represents the network of the last 100 tweets using the hashtag #aacrao10. The data was captured and the network was constructed using R. Thanks to Drew Conway for showing me how to do this.

(Click to enlarge)



Labeled dots indicate users that used the specified hashtag while unlabeled dots indicate 'friends' of users that used the specified hashtag.The code also allows you to compute important network metrics such as measures of centrality that are helpful in key actor analysis.

Statistics Definitions

SAMPLE POINT- each individual outcome of an experiment

SAMPLE SPACE-the collection of all possible sample points in an experiment

EVENT-a collection of sample points

PRIOR PROBABILITY-initial estimate of the probability of an event

RANDOM VARIABLE-a numerical description of the outcome of an experiment

EXPECTED VALUE OF A RANDOM VARIABLE-a measure of the average value of a random variable

VARIANCE OF A RANDOM VARIABLE-a measure of the dispersion of a random variable

PARAMETER-a numerical measure from a population

STATISTIC- a numerical measure from a sample

SAMPLING DISTRIBUTION-a probability distribution for all possible values of a statistic

PROPERTIES OF ESTIMATORS:

Unbiased: Property of a statistic where the expected value of a statistic/estimator = the parameter being estimated
Consistency: as n increases the probability that the value of a statistic/estimator gets closer to the parameter being estimated increases

CENTRAL LIMIT THEOREM-the sampling distribution of the sample mean can be approximated by a normal probability distribution as the sample size becomes large. This holds even if the population from which the sample mean comes from is not normal, or unknown.

ASYMPTOTIC DISTRIBUTION- an estimator’s sampling distribution in large samples

ASYMPTOTIC PROPERTIES OF AN ESTIMATOR- properties of an estimator in large samples