Lab4
Priya Shaji
3/2/2019
A) Foundations for statistical inference - Sampling distributions
Lab
The data
Download the file and Load the data
#library(DATA606)
#download.file("http://www.openintro.org/stat/data/ames.RData", destfile = "ames.RData")
load("/Users/priyashaji/Documents/cuny msds/Spring'19/data 606/Labs/Lab4/ames.RData")Renaming columns to be used
area <- ames$Gr.Liv.Area
price <- ames$SalePriceLet’s look at the distribution of area in our population of home sales by calculating a few summary statistics and making a histogram.
summary(area)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 334 1126 1442 1500 1743 5642hist(area)
mean(area)## [1] 1499.69The unknown sampling distribution
If we were interested in estimating the mean living area in Ames based on a sample, we can use the following command to survey the population.
samp1 <- sample(area, 50)
summary(samp1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 816 1040 1440 1432 1675 3500hist(samp1)
If we’re interested in estimating the average living area in homes in Ames using the sample, our best single guess is the sample mean.
mean(samp1)## [1] 1431.86Generate 5000 samples and compute the sample mean of each.
sample_means50 <- rep(NA, 5000)
for(i in 1:5000){
samp <- sample(area, 50)
sample_means50[i] <- mean(samp)
}
hist(sample_means50)
If you would like to adjust the bin width of your histogram to show a little more detail, you can do so by changing the breaks argument.
hist(sample_means50, breaks = 25)
Interlude: The for loop
Calculating the sample mean without the for loop
With the for loop, these thousands of lines of code are compressed into a handful of lines. We’ve added one extra line to the code below, which prints the variable i during each iteration of the for loop. Run this code.
Sample size and the sampling distribution
Mechanics aside, let’s return to the reason we used a for loop: to compute a sampling distribution, specifically, this one.
hist(sample_means50)
To get a sense of the effect that sample size has on our distribution, let’s build up two more sampling distributions: one based on a sample size of 10 and another based on a sample size of 100.
sample_means10 <- rep(NA, 5000)
sample_means100 <- rep(NA, 5000)
for(i in 1:5000){
samp <- sample(area, 10)
sample_means10[i] <- mean(samp)
samp <- sample(area, 100)
sample_means100[i] <- mean(samp)
}To see the effect that different sample sizes have on the sampling distribution, plot the three distributions on top of one another.
par(mfrow = c(3, 1))
xlimits <- range(sample_means10)
hist(sample_means10, breaks = 20, xlim = xlimits)
hist(sample_means50, breaks = 20, xlim = xlimits)
hist(sample_means100, breaks = 20, xlim = xlimits)
Exercises
Exercise 1
Describe this population distribution.
Answer 1)
hist(area)
The distribution is left skewed and unimodal
Exercise 2
Describe the distribution of this sample. How does it compare to the distribution of the population?
samp1 <- sample(area, 50)
summary(samp1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 848 1221 1474 1500 1676 2550mean(samp1)## [1] 1499.8hist(samp1)
The distribution of the sample seems to be unimoadl and left skewed. It looks similar to the distribution of the population.
Exercise 3
Take a second sample, also of size 50, and call it samp2. How does the mean of samp2 compare with the mean of samp1? Suppose we took two more samples, one of size 100 and one of size 1000. Which would you think would provide a more accurate estimate of the population mean?
samp2 <- sample(area, 50)
summary(samp2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 660 1008 1339 1409 1694 2646mean(samp2)## [1] 1409.28hist(samp2)
Mean of sample 1: 1537.16
Mean of sample 2: 1616.56
Sample size : 100
samp3 <- sample(area, 100)
summary(samp3)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 672 1051 1464 1469 1729 3228mean(samp3)## [1] 1469.39hist(samp3)
Sample size: 1000
samp4 <- sample(area, 1000)
summary(samp4)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 334 1120 1454 1506 1759 5642mean(samp4)## [1] 1505.757hist(samp4)
The mean of the samp2 will not be the same as the mean of samp1. It is difficult to determine if samp2 will consistantly have a higher or lower mean because the samples are randomly selected.
The larger the sample mean, the more accurate the estimate of the population mean will be. As such, the sample size of 1,000 will be the closest estimate to the original population.
Exercise 4
How many elements are there in sample_means50? Describe the sampling distribution, and be sure to specifically note its center. Would you expect the distribution to change if we instead collected 50,000 sample means?
sample_means50 = rep(NA, 5000)
for(i in 1:5000){
samp = sample(area, 50)
sample_means50[i] = mean(samp)
}
hist(sample_means50, breaks = 25, main = "Area Sample_Means_50")
The sampling distribution is normal and unimodal with a center of approximately 1500. If sample means is 50,000, the more accurate the estimate of the population mean will be. And sample size of 1,000 will be the closest estimate to the original population.
Exercise 5
To make sure you understand what you’ve done in this loop, try running a smaller version. Initialize a vector of 100 zeros called sample_means_small. Run a loop that takes a sample of size 50 from area and stores the sample mean in sample_means_small, but only iterate from 1 to 100. Print the output to your screen (type sample_means_small into the console and press enter). How many elements are there in this object called sample_means_small? What does each element represent?
sample_means_small = rep(NA, 100)
for(i in 1:100){
samp = sample(area, 50)
sample_means_small[i] = mean(samp)
}
sample_means_small## [1] 1552.98 1536.86 1548.36 1472.72 1489.04 1545.72 1496.46 1501.22
## [9] 1420.24 1436.16 1448.12 1467.76 1549.64 1447.10 1509.54 1482.94
## [17] 1490.00 1658.96 1463.72 1592.54 1422.22 1515.28 1460.80 1470.80
## [25] 1515.18 1525.96 1486.72 1488.50 1467.10 1499.50 1461.34 1446.32
## [33] 1487.78 1426.98 1401.00 1336.68 1433.40 1615.28 1484.74 1542.08
## [41] 1355.94 1457.84 1447.58 1566.34 1519.08 1473.00 1485.54 1458.34
## [49] 1656.36 1411.96 1451.48 1459.12 1647.18 1463.72 1472.02 1528.42
## [57] 1526.70 1471.30 1575.02 1473.56 1497.06 1526.20 1629.82 1453.46
## [65] 1598.66 1482.92 1516.96 1478.76 1423.36 1527.28 1557.84 1485.68
## [73] 1462.12 1529.88 1462.98 1499.80 1452.98 1670.62 1524.60 1454.72
## [81] 1530.62 1373.22 1399.50 1333.14 1482.76 1479.60 1548.70 1479.54
## [89] 1499.94 1323.82 1447.36 1512.74 1523.24 1406.16 1339.02 1476.88
## [97] 1567.16 1451.86 1550.30 1524.80hist(sample_means_small, breaks = 25, main = "Area Sample_means_small")
There are 100 elements in sample_means_small object.
Each element represents the mean of a sample of size 50, drawn for the population area.
Exercise 6
When the sample size is larger, what happens to the center? What about the spread?
When the sample size is larger, the spread decreases and the center grows closer to the population area’s center.
On your Own
Ques_1
Take a random sample of size 50 from price. Using this sample, what is your best point estimate of the population mean?
price <- ames$SalePrice
samp1 <- sample(price, 50)
summary(samp1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 76500 121250 152750 175073 199250 395039mean(samp1)## [1] 175073.1hist(samp1)
The mean is the best point estimate for the population mean.
Ques_2
Since you have access to the population, simulate the sampling distribution for x¯price by taking 5000 samples from the population of size 50 and computing 5000 sample means. Store these means in a vector called sample_means50. Plot the data, then describe the shape of this sampling distribution. Based on this sampling distribution, what would you guess the mean home price of the population to be? Finally, calculate and report the population mean.
sample_means50 = rep(NA, 5000)
for(i in 1:5000){
samp = sample(price, 50)
sample_means50[i] = mean(samp)
}
summary(sample_means50)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 146065 172984 180509 180963 188507 223795hist(sample_means50, breaks = 25, main = "Price Sample_means50")
Shape of the sampling distribution is normal and unimodal. The mean home price of the population seems to be 181121.
Calculating mean for population:
summary(price)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12789 129500 160000 180796 213500 755000The actual mean of the price is 180796. My estimate was higher than the actual mean but much closer to the actual value.
Ques_3
Change your sample size from 50 to 150. Repeat process above but store these means in a vector called “sample_means150”. Describe the shape of this sampling distribution. Compare it to the sampling distribution for a sample size of 50. Based on this sampling distribution, what would you guess to be the mean sale price of homes in Ames?
sample_means150 = rep(NA, 5000)
for(i in 1:5000){
samp = sample(price, 150)
sample_means150[i] = mean(samp)
}
summary(sample_means150)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 157275 176480 180616 180816 185077 209761hist(sample_means150, breaks = 25, main = "Price Sample_means50")
The sample distribution is normal and unimodal.
# Comparisons between means50 and means150
par(mfrow= c(2,1))
hist(sample_means50, breaks = 50, main = "Price Sample_Means_50")
hist(sample_means150, breaks = 50, main = "Price Sample_Means_150")
summary(sample_means50)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 146065 172984 180509 180963 188507 223795summary(sample_means150)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 157275 176480 180616 180816 185077 209761The sampling distribution for sample size 50 is more spread as compared to the sample size of 150.
The mean sale price seems to be 180725, which is approx the actual mean sale price of homes in ames.
Ques_4
Of the sampling distributions from 2 and 3, which has a smaller spread? If we’re concerned with making estimates that are more often close to the true value, would we prefer a distribution with a large or small spread?
par(mfrow= c(2,1))
hist(sample_means50, breaks = 50, main = "Price Sample_Means_50")
hist(sample_means150, breaks = 50, main = "Price Sample_Means_150")
Of the two sampling distributions (sample_means50 and sample_means150), the second (sample_means150) has the smaller spread.
It is preferable to have a distribution with a smaller spread when trying to make estimates that are close to the true value, because there is less uncertainty of where the true estimates lie.
B) Foundations for statistical inference - Confidence intervals
Lab
The data
Download the file and Load the data
#library(DATA606)
#download.file("http://www.openintro.org/stat/data/ames.RData", destfile = "ames.RData")
load("/Users/priyashaji/Documents/cuny msds/Spring'19/data 606/Labs/Lab4/ames.RData")Let’s start with a random sample size of 60
population <- ames$Gr.Liv.Area
samp <- sample(population, 60)Confidence intervals
One of the most common ways to describe the typical or central value of a distribution is to use the mean. In this case we can calculate the mean of the sample using,
sample_mean <- mean(samp)We can calculate a 95% confidence interval for a sample mean by adding and subtracting 1.96 standard errors to the point estimate (See Section 4.2.3 if you are unfamiliar with this formula).
se <- sd(samp) / sqrt(60)
lower <- sample_mean - 1.96 * se
upper <- sample_mean + 1.96 * se
c(lower, upper)## [1] 1425.031 1672.802Confidence levels
In this case we have the luxury of knowing the true population mean since we have data on the entire population. This value can be calculated using the following command:
mean(population)## [1] 1499.69we need to first create empty vectors where we can save the means and standard deviations that will be calculated from each sample. And while we’re at it, let’s also store the desired sample size as n.
samp_mean <- rep(NA, 50)
samp_sd <- rep(NA, 50)
n <- 60Now we’re ready for the loop where we calculate the means and standard deviations of 50 random samples.
for(i in 1:50){
samp <- sample(population, n) # obtain a sample of size n = 60 from the population
samp_mean[i] <- mean(samp) # save sample mean in ith element of samp_mean
samp_sd[i] <- sd(samp) # save sample sd in ith element of samp_sd
}Lastly, we construct the confidence intervals.
lower_vector <- samp_mean - 1.96 * samp_sd / sqrt(n)
upper_vector <- samp_mean + 1.96 * samp_sd / sqrt(n)Lower bounds of these 50 confidence intervals are stored in lower_vector, and the upper bounds are in upper_vector. Let’s view the first interval.
c(lower_vector[1], upper_vector[1])## [1] 1394.019 1613.314Exercises
Exercise 1
Describe the distribution of your sample. What would you say is the “typical” size within your sample? Also state precisely what you interpreted “typical” to mean.
hist(samp,breaks = 50)
The population unimodal and left skewed.
The typical size within a sample ranges within 1300 to 1450
By “typical size” I analyzed the size which has highest number of frequency.
Exercise 2
Would you expect another student’s distribution to be identical to yours? Would you expect it to be similar? Why or why not?
Another student’s distribution will bw different,as the sample of population might vary for different students and also it is relatively small sample size.
The sample size is 60 out of 5000 cases.The distributions can be similar on most rare occasions.
Exercise 3
For the confidence interval to be valid, the sample mean must be normally distributed and have standard error s/√‾n. What conditions must be met for this to be true?
The following condtions should be met:
The observations should be independent.
The sample size should be more than 30.
The sample data should not be strongly skewed.
Exercise 4
What does “95% confidence” mean?
It means that we are 95 % sure about the true population mean lying between interval (-2 Sd - +2 Sd) when ploting sampling distribution. See the range below:
lower## [1] 1425.031upper## [1] 1672.802Exercise 5
Does your confidence interval capture the true average size of houses in Ames? If you are working on this lab in a classroom, does your neighbor’s interval capture this value?
summary(population)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 334 1126 1442 1500 1743 5642Yes, if you see the mean of the population above, it does fall under lower and upper level of the sample calculated using 95% confidence.
Exercise 6
Each student in your class should have gotten a slightly different confidence interval. What proportion of those intervals would you expect to capture the true population mean? Why? If you are working in this lab in a classroom, collect data on the intervals created by other students in the class and calculate the proportion of intervals that capture the true population mean.
95% of the students will capture true population because of the fact that everyone of them would have used the following (Sample Mean + (+/- 1.96 * SE))
On your own
Ques_1
Using the following function (which was downloaded with the data set), plot all intervals. What proportion of your confidence intervals include the true population mean? Is this proportion exactly equal to the confidence level? If not, explain why.
plot_ci(lower_vector, upper_vector, mean(population))
In the first run I got 46/50 i.e. 92% which is lower than 95%. The 95% confidence level does not mean that every run would have poopulation mean lying between the sample mean interval. There would be cases where the population mean would be lying outside interval (Sample Mean + (+/- 1.96 * SE)).
Ques_2
Pick a confidence level of your choosing, provided it is not 95%. What is the appropriate critical value?
We can choose the confidence level of 99% by following the given formula: (Sample Mean + (+/- 2.58 * SE)).
Ques_3
Calculate 50 confidence intervals at the confidence level you chose in the previous question. You do not need to obtain new samples, simply calculate new intervals based on the sample means and standard deviations you have already collected. Using the plot_ci function, plot all intervals and calculate the proportion of intervals that include the true population mean. How does this percentage compare to the confidence level selected for the intervals?
lower_vector <- samp_mean - 2.58 * samp_sd / sqrt(n)
upper_vector <- samp_mean + 2.58 * samp_sd / sqrt(n)
plot_ci(lower_vector, upper_vector, mean(population))