Data Visualisation 2019 - Assignment 4

List of R colors:

http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf

The normal distribution

• The normal distribution is the graph of the function \[ f(x)=\frac{1}{\sqrt{2\pi \sigma}}e^{-\frac{(x-\mu)^2}{2\sigma}}\] where \(\mu\) (pronounced as “mu”) is the mean and \(\sigma\) (pronounced as “sigma”") is the standard deviaiton of the distribution.

Example 1

• The standard normal distribution is the one where \(\mu=0\) and \(\sigma=1\), and is shown in the graph below

X<-seq(-5,5,len=201) #Create the x-axis between -5 and +5
P<-dnorm(X,mean=0,sd=1) # Evaluate the Normal Distribution function over the x-axis
plot(X,P,type='l',lwd=2,col='steelblue') # Plot the Probability along th x-axis

• The lines of this cell have the following functions:

1. X<-seq(-5,5,len=200) This creates a sequence of \(x\)-values, which will form the horizontal axis. These values start at -5, and stop at +5, and the sequence has 200 points altogether.

X  

  [1] -5.00 -4.95 -4.90 -4.85 -4.80 -4.75 -4.70 -4.65 -4.60 -4.55 -4.50 -4.45 -4.40 -4.35 -4.30 -4.25 -4.20 -4.15 -4.10 -4.05 -4.00 -3.95 -3.90 -3.85 -3.80 -3.75 -3.70
 [28] -3.65 -3.60 -3.55 -3.50 -3.45 -3.40 -3.35 -3.30 -3.25 -3.20 -3.15 -3.10 -3.05 -3.00 -2.95 -2.90 -2.85 -2.80 -2.75 -2.70 -2.65 -2.60 -2.55 -2.50 -2.45 -2.40 -2.35
 [55] -2.30 -2.25 -2.20 -2.15 -2.10 -2.05 -2.00 -1.95 -1.90 -1.85 -1.80 -1.75 -1.70 -1.65 -1.60 -1.55 -1.50 -1.45 -1.40 -1.35 -1.30 -1.25 -1.20 -1.15 -1.10 -1.05 -1.00
 [82] -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05  0.00  0.05  0.10  0.15  0.20  0.25  0.30  0.35
[109]  0.40  0.45  0.50  0.55  0.60  0.65  0.70  0.75  0.80  0.85  0.90  0.95  1.00  1.05  1.10  1.15  1.20  1.25  1.30  1.35  1.40  1.45  1.50  1.55  1.60  1.65  1.70
[136]  1.75  1.80  1.85  1.90  1.95  2.00  2.05  2.10  2.15  2.20  2.25  2.30  2.35  2.40  2.45  2.50  2.55  2.60  2.65  2.70  2.75  2.80  2.85  2.90  2.95  3.00  3.05
[163]  3.10  3.15  3.20  3.25  3.30  3.35  3.40  3.45  3.50  3.55  3.60  3.65  3.70  3.75  3.80  3.85  3.90  3.95  4.00  4.05  4.10  4.15  4.20  4.25  4.30  4.35  4.40
[190]  4.45  4.50  4.55  4.60  4.65  4.70  4.75  4.80  4.85  4.90  4.95  5.00

2. P<-dnorm(X,mean=0,sd=1) This evaluates the normal distribution with mean=0 and standard deviation = 1. With mean=0 this means the peak of the curve is at \(x=0\). The \(sd=1\) command determines how spread out the graph is, the larger sd the more spread out the graph.

3. The last line instruct R to plot X on the horizontal axis and P on the vertical axis. The argument type=‘l’ instructs R to plot the curve as a line, lwd=2 indicates a linewidth of 2 and col=‘steelblue’ indicates the colour to plot the curve.

Exercise 1

Plot the normal disribution with \(\mu=4\) and \(\sigma=2.15\). Note that you will have to re-define an appropriate X-interval for this plot also. A convient way to do this so that your plot looks symmetric is to set \[\mathbf{X=\text{seq}(\mu-5*\sigma,\mu+5*\sigma)}\] and do not forget to include * when you want to denote multiplication.

Interpretation of the normal curve

• When \(x\) is normally distributed with mean \(\mu\) and standard deviation \(\sigma\), then the proability that \(x\) lies between the values \(x_1\) and \(x_2\) is the area under the normal distribution curve, between these two values.

• An improtant feature of every normal distribution curve is the the total area beneath the curve is always equal to 1.

Example 2

Illustrate the probability for \(x\) in the intervavl \(-2\leq x\leq 1\) in the normal distribution with \(\mu=0\) and \(\sigma=1\). Calculate this proability.

# Part 1. Plot the normal curve again - See Example 1
X<-seq(-5,5,len=201)
P<-dnorm(X,mean=0,sd=1)
plot(X,P,type='l',lwd=2,col='steelblue')

# Part 2. Shade the required area.
X1=seq(-2,1,length=201) # The interval -5<x<5
P1=dnorm(X1,mean=0,sd=1)
polygon(c(-2,X1,1),c(0,P1,0),col=adjustcolor("cornflowerblue",alpha.f=0.3),border=NA)

• The code in the cell above has the following meaning:

1. Part 1 This code has alredy been explained in Example 1

2. Part 2 We create a new x1 axis between -5 and -1. We create the corresponding y-values using the dnorm() funciton applied to the intrval x1. polygon(c(-5,x_1,-1),c(0,y1,0)) indicates that the area between -5__ and -1 on the x-axis and between 0 and the corresponging values of y1 should be filled in.

3. col=adjustcolor(“cornflowerblue”,alpha.f=0.3) indicates the fill color and the transparencey alpha of the filling. The smaller the value of alpha the more transparent the filling, where alpha should always be between 0 and 1.

Probabilities

• The probability that \(x\) lies between \(-\infty\) and \(-1\) is found using the function dnorm()

dnorm(-1,mean=0,sd=1)

[1] 0.2419707

• The probability that \(x\) lies between \(-\infty\) and \(2\) is

dnorm(2,mean=0,sd=1)

[1] 0.05399097

• The the probability that x lies between -1 and 2 is

dnorm(2,mean=0,sd=1)-dnorm(-1,mean=0,sd=1)

[1] -0.1879798

• This correspond to the shaded area shown above.

Confidence Intervals

Given a population, a Confidence Interval (CI) is a means of estimating a population statistic from the corresponding sample staistic with a certain level of confidence. For example, we may take survey data on a collection of people to estimate the average mothly salary of the people in the survey. Then based on this sample average we can estimate a range for the average slalary of all people in the population, to a certain confidence level. The range of possible values for the population average is called a confidence interval.

Example 3:

• Suppose we wish to determine the mean salary of the population of a country with a certain degree of confidence.

• In practice it would be difficult, if not impossible, to collect data on the salary of everybody from the population.

• Instead, we may collect salary data from numerous samples of people from the population, and then get the mean salary for each of these samples.

• For each sample, the sample mean would be spread with a normal distribution, whose standard deviation and mean would be close to the mean of the overall population.

• We let \(\bar{x}\) denote the mean salary from one of these samples, and the standard deviation of the sample be \(\sigma\) (this measures how much the means are spread about the true mean).

• We can now say that the true population mean lies in the interval \[\bar{x}\pm z_{\frac{\alpha}{2}}\frac{\sigma}{\sqrt{n}},\] where \(n\) is the size of the sample for which the mean \(\bar{x}\) was calculated.

• The parameter \(\alpha\) is called the confidence parameter associated with the confidence level, and the number \(z_{\frac{\alpha}{2}}\) is the z-score associated with the parameter \(\frac{\alpha}{2}\).

• From a sample of 100 people, we find that the lowest salary from the sample is EUR 15000 and the highest is EUR 55000, with a mean salary of EUR 35000 and a standard deviation of EUR 2500.

• We can plot the normal distribution for this data as follows:

  1. We create a sample size between 15000 and 55000 (the sample range), with step size 50. We call this sample data Sample1

  2. We use the dnorm() function to indicate that the data in Sample1 is normally distributed about the mean where the mean is 35000 and the standard deviation is 2500.

  3. We plot the distribution using the usual plot() function, with Sample1 on the x-axis and Data1 on the y-axis.

Sample1<-seq(15000,55000, by=50)
Data1 <- dnorm(Sample1,mean=35000, sd=2500)
plot(Sample1,Data1,ylab="Rel. Freq.",xlab="Salaries (EUR)",type="l",lwd=2,col="red")

A 90% Confidence Interval

• Suppose we wish to determine the average salary value at a confidence level of 90%.

• This would correspond to the salary interval which gives 90% of the area beneath the normal distribution curve above.

• This number can be calculated using a normal distribution table, and finding the \(z_{\alpha}\)-score corresponding to \(\alpha=1-0.90=0.1\).

# Part 1. Plot the normal curve again - See Example 1
Z0<-seq(-5,5,len=200)
H0<-dnorm(Z0,mean=0,sd=1)
plot(Z0,H0,type='l',lwd=2,col='steelblue',main='The 90% confidence interval')

# Part 2. Shade the required area.
x1=seq(-qnorm(0.95),qnorm(0.95),length=200)
y1=dnorm(x1,mean=0,sd=1)
polygon(c(-qnorm(0.95),x1,qnorm(0.95)),c(0,y1,0),col=adjustcolor("cornflowerblue",alpha.f=0.3),border=NA)

• When we look for a confidence interval of 90% we are looking for the values of \(z\) such that 90% of the area beneath the curve, between \(-z_{\frac{\alpha}{2}}\) and \(z_{\frac{\alpha}{2}}\).

• It also means we are looking for the \(z\)-values which will exclude 5% of the area in each tail, that is we are looking for the \(z\)-values for each of the dashed line in the plot below:

Standard<-seq(-4,4,by=0.001)
StandardData<-dnorm(Standard,mean=0,sd=1)
plot(Standard,StandardData,ylab="Rel. Freq.", xlab="z", type="l", lwd=2,col="red")
segments(-qnorm(0.9),0,-qnorm(0.9),dnorm(-qnorm(0.9),mean=0,sd=1),lwd=2.5,lty=2)

segments(qnorm(0.9),0,qnorm(0.9),dnorm(qnorm(0.9),mean=0,sd=1),lwd=2.5,lty=2)

• 90% of the area beneath the curve is contained between the two dashed lines.

• 5% of the area is in each tail to the left and right.

• It also means that the area beneath the curve containing the left tail and the middle portion is 95%

• To find the \(z\)-value corresponding to this 95 % in R we can use the qnorm() function as follows:

z90=qnorm(0.95)
z90

[1] 1.644854

• In the plot above this means that 95% of the are beneath the curve is between \(z=-\infty\) and \(1.64485\).

• It also means that 90% of the area lies, from \(z=-1.644854\) to \(z=-1.644854\).

• For the purposes of our salary survey above, it also means \[z_{\frac{\alpha}{2}}=1.644845\]

• In terms of our survey we are also 90% certain that the mean population salry lies in the range

lower_bound=35000-z90*2500/sqrt(100)
upper_bound=35000+z90*2500/sqrt(100)
lower_bound

[1] 34588.79

upper_bound

[1] 35411.21

We are 90% confident that the avareage salary of the overall population lies in the range €34,588.79 to €35,411.21.

In terms of our normal distribution plot this means, we are 90% certain the mean population salary lies between the two dashed blue lines.

Sample1<-seq(15000,55000, by=50)
Data1 <- dnorm(Sample1,mean=35000, sd=2500)
plot(Sample1,Data1,ylab="Rel. Freq.",xlab="Salaries (EUR)",type="l",lwd=2,col="red")
abline(v=lower_bound,lwd=2,lty=2,col="blue")

abline(v=upper_bound,lwd=2,lty=2,col="blue")

Exercise 2

A company wishes to determine the cost associated with providing health insurance for all of its employees. The company asked a sample of 35 employees how much they currently pay annually for health insurance. The data obtained revealed a sample mean premium of EUR635 with a standard deviation of EUR67. Using this data answer the following

1. Find the 90%, 95%, 96% and 99% confidence intervals for the average amount spent by all customers.
2. Plot a normal distribution for this data set showing the confidence intervals obtained.

Exercise 3

The mean number of emails received by employees at a company was measured on a ample of 58 employees. The sample mean was found to be 21 emails per day with a standard deviation of 9. Using this data answer the following

1. Find the 90%, 94%, 96% and 99% confidence intervals for the average amount spent by all customers.
2. Plot a normal distribution for this data set showing the confidence intervals obtained.

Samples, Means and Standard Deviations

The Mean

Given a data set we can calculate the mean and the standard deviation of the data-set using the functions mean() and sd().

Given a data set of size \(n\), with data values \(x_1,x_2,\ldots,x_n\), the mean of the data set is given by \[ \bar{x}=\frac{\displaystyle{\sum_{i=1}^{n}}x_{i}}{N}. \] This can be evaluated in R using the function mean():

Example 2

Suppose the data from a survey is given by the data set \[ 1.3,-4.2,1.55,6.78,4.5,3.21, 5.89,7.34,-6.88,-7.9 \] Create a data structure in R to represent this data set. Use the mean() function to evaluate the mean of the data set.

dataset2<-c(1.3,-4.2,1.55,6.78,4.5,3.21, 5.89,7.34,-6.88,-7.9)
mean(dataset2)

[1] 1.159

• Hence the mean is \[\bar{x}=1.159\]

The Standard Deviation

Given a data set we can calculate the mean and the standard deviation of the data-set using the functions mean() and sd().

Given a data set of size \(n\), with data values \(x_1,x_2,\ldots,x_n\), the mean of the data set is given by \[ s=\frac{\displaystyle{\sum_{i=1}^{n}}(x_{i}-\bar{x})^2}{N-1}. \] This can be evaluated in R using the function mean():

Example 4

Evaluate the standard deviation of the data set in Example 2 using the sd() function.

sd(dataset2)

[1] 5.608865

• Hence the standard deviation is \[ s= 5.608865 \]

Example 5

The inflation rate in Ireland between the years 1992 and 2016 is given in the data file IrelandInflation.csv(Source: World Bank). We will import this into R using the usual command read.csv(file.choose()) as a data structure. Using this data structure we want to answer the following:

1. Find the mean inflation rate between the years 1992 and 2016
2. Find the standard deviation of this data set
3. Find the number of data points in this data set
4. Using these three values, find the range of possible values for the average inflation rate in Ireland over all years with 95% confidence
5. Plot the inflation data as a normal distribution.

Importing the Data

InflationData<-read.csv("IrelandInflation.csv")
InflationData

Creating the Data Structure

Inflation<-InflationData$Inflation.Rate..
Inflation

 [1]  3.119790e+00  1.408776e+00  2.345707e+00  2.514464e+00  1.693076e+00  1.437211e+00  2.426878e+00  1.640271e+00  5.564830e+00  4.872355e+00  4.651952e+00
[12]  3.479883e+00  2.194873e+00  2.431541e+00  3.938895e+00  4.879925e+00  4.053506e+00 -4.479938e+00 -9.461664e-01  2.578870e+00  1.692785e+00  5.026782e-01
[23]  1.967858e-01 -2.945990e-01  3.120000e-13

Part 1 - The Mean Inflation

xbar=mean(Inflation)
xbar

[1] 2.076174

• The mean inflation rate between 1992 and 2016 in Ireland was \(\bar{x}=2.076174%\)

Part 2 - The Standard Deviation

s=sd(Inflation)
s

[1] 2.189926

• The standard deviation in the rate of inflation between 1992 and 2016 is \(\sigma=2.189926\)

Part - The Number of Data Points

We use the length() function to determine the number of data-points in Inflation

n=length(Inflation)
n

[1] 25

• There are \(n=25\) data points in the data set Inflation

Part 4 - The 90% Confidence Interval

• The confidence parameter \(\alpha\) is given by \[\alpha = 1-0.95=0.05\Rightarrow \frac{\alpha}{2}=0.025\]

• We now wish to find the \(z_{\frac{\alpha}{2}}\) which gives an area of \(0.95+0.025=0.975\) under the curve of the normal distribution. We find this \(z\)-score using qnorm()

z0=qnorm(0.975)
z0

[1] 1.959964

• For the 95% confidence interval, we have \(z_{\frac{\alpha}{2}}=1.959964\). Hence the upper and lower bound for the 95% confidence interval for the data set are

xlower=xbar-z0*s/sqrt(n)
xupper=xbar+z0*s/sqrt(n)

xlower

[1] 1.217739

xupper

[1] 2.934609

We are 95% certain that the mean inflation rate for Ireland lies in the interval 1.217739% to 2.934609%.

Part 5 - The Normal Distribution Plot

We now plot StandardInflation to give the normal distribution plot:

1. First we create an interval appropriate for creating a normal distribution form the data-set Inflation: For this plot, we will make the plot range start at \(\bar{x}-5\sigma\) and end at \(\bar{x}+5\sigma\), with a step-size 0.0001:

Interval1=seq(xbar-5*s,xbar+5*s,by=0.0001)

2. Next we create a normal distribution from this data with mean and standard deviation evaluated earlier:

NormalDistInflation=dnorm(Interval1, mean=xbar, sd=s)

3 Now we plot NormalDistInflation on the y-axis and Interval1 on the x-axis:

plot(Interval1,NormalDistInflation,ylab="Rel. Freq.",xlab="Inflation Rate %",type="l",lwd=2,col="red")
abline(v=xlower,lwd=2,lty=2,col="blue")

abline(v=xupper,lwd=2,lty=2,col="blue")

Exercise 4:

• The data file IrelandEmissions.csv contains the CO2 emissions, in metric tons per capital, of Ireland between 1992 and 2016. (Source: World Bank)

• Using the data in this file, create a data structure to capture the numerical data over those 25 years.

• Using this data structure answer the following:

1. Find the mean inflation rate between the years 1992 and 2016
2. Find the standard deviation of this data set
3. Find the number of data points in this data set
4. Using these three values, find the 90% and 99% confidence intervals for the average CO2 emissions of Ireland over all years.
5. Plot the inflation data as a normal distribution, and indicate the confidence intervals on these plots.

Interpretation of Confidence Intervals

• We record some data associated with a population, given in the data set X below:

X=c(0.34,0.77,0.98,1.03,1.45,1.87,1.99,2.11,2.23,2.24,2.87,2.89,2.97,3.04,3.45,3.56,3.77,3.79,3.81,3.97,4.12,4.34,4.88,5.01,5.16)

• Obviously, the population mean is easily available to us in this case:

mean(X)

[1] 2.9056

Sampling the population

• In many cases of practical interest, this mean is not available to us, and the best we can do is form a confidence interval for this mean.

• To form these confidence intervals, we take samples of the population data.

• This is implemented in R as follows:

n=4 # sample size
N=10 # number of samples
Sample<-replicate(N,sample(X,n))

• We have taken 10 Samples with each sample having 4 data values

Sample

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 3.45 2.24 1.99 1.03 0.77 0.98 3.56 3.97 5.01  1.87
[2,] 3.77 5.16 4.12 0.34 2.24 3.45 5.16 3.56 2.89  3.97
[3,] 2.23 2.89 2.23 2.23 2.97 2.97 1.45 0.34 3.45  1.03
[4,] 1.99 1.03 2.87 5.16 3.97 3.04 2.97 4.88 3.77  4.88

• Now we wish to take the mean of each sample, which we do as follows:

M <- numeric(length = length(1:N))
for (i in 1:N) {
  M[i] <- mean(Sample[,i])
}
M

 [1] 2.8600 2.8300 2.8025 2.1900 2.4875 2.6100 3.2850 3.1875 3.7800 2.9375

• Clearly, the mean of each sample does not match the mean of the population.

Confidence Intervals from the Samples

• We now form the 60% Confidence Interval associated with each sample mean, and to do so we need the z score corresponding to the 80th percentile:

z<-qnorm(0.80,mean=0,sd=1)
z

[1] 0.8416212

• Next we form the standard error associated with the population and the samples:

1. The standard deviation of the population is

sigma<-sqrt(sum((X-mean(X))^2)/length(X))
sigma

[1] 1.343143

We do this because the function sd() uses the formula for sample standard deviation.

2. The standard error is

SE<-sigma/sqrt(n)
SE

[1] 0.6715714

i.e. the population standard deviation divided by the sample size.

3. The lower boundaries of the confidence intervals are now:

ML<-M-z*SE
ML

 [1] 2.294791 2.264791 2.237291 1.624791 1.922291 2.044791 2.719791 2.622291 3.214791 2.372291

The upper ends of the confidence intervals are:

MU<-M+z*SE
MU

 [1] 3.425209 3.395209 3.367709 2.755209 3.052709 3.175209 3.850209 3.752709 4.345209 3.502709

• The 60% confidence intervals following from these 10 samples are:

paste(ML,MU,sep=' to ')

 [1] "2.29479124068733 to 3.42520875931267" "2.26479124068733 to 3.39520875931267" "2.23729124068733 to 3.36770875931267" "1.62479124068733 to 2.75520875931267"
 [5] "1.92229124068733 to 3.05270875931267" "2.04479124068733 to 3.17520875931267" "2.71979124068733 to 3.85020875931267" "2.62229124068733 to 3.75270875931267"
 [9] "3.21479124068733 to 4.34520875931267" "2.37229124068733 to 3.50270875931267"

Recall: The actual mean of the pupulation is

mean(X)

[1] 2.9056

• We deliberitely chose 60% confidence intervals to highlight the point that not all confidence intervals contain the true population mean.

• This gives us the correct interpretation for the 60% confidence interval:

  • A confidence iterval either does or does not contain the true population mean.

  • As we can see, if we take enough samples, then some of these confidence intervals will contain the mean and some will not.

  • If we take a large enough number of samples, then on average, 60% of the confidence intervals deduced from those samples, at confidence level 60%, will contain the true mean.

Exercise 5:

1. Using the same data set X, take N=50 samples of sample size n=5.

2. Find the mean of each sample.

3. Find the 90% confidence interval associated with these sample averages.

4. Do all these confidence intervals contain the true population mean?