## CIs for mu and t-tests about mu:

###################################################################################################
###################################################################################################
###################################################################################################
#
# Copy and paste this section into R:
#
#
t.test.stats <- function(x.bar, s, n, mu=0, alternative="two.sided"){
  t.stat <- (x.bar-mu)/(s/sqrt(n))
  if (alternative=="less") P.value <- pt(t.stat,df=n-1,lower=T)
  if (alternative=="greater") P.value <- pt(t.stat,df=n-1,lower=F)
  if (alternative=="two.sided") P.value <- 2*(min(pt(t.stat,df=n-1,lower=T), pt(t.stat,df=n-1,lower=F)))
  my.out <- list(t.stat, P.value, alternative)
  names(my.out) <- c("test statistic value", "P-value", "alternative")
  return(my.out)
}
#
t.CI.stats <- function(x.bar, s, n, conf.level=0.95){
  alpha<- 1-conf.level
  lower <- x.bar - qt(alpha/2,df=n-1,lower=F)*(s/sqrt(n))
  upper <- x.bar + qt(alpha/2,df=n-1,lower=F)*(s/sqrt(n))
  my.out <- list(lower, upper, conf.level)
  names(my.out) <- c("lower limit", "upper limit", "confidence level")
  return(my.out)
}
#
#
###################################################################################################
###################################################################################################
###################################################################################################

## Example 1:

########## If you are given the entire data set:

my.datafile <- tempfile()
cat(file=my.datafile, "
1 Austria 8.6
2 Belgium 10.3
3 Czechrepublic 11.3
4 Denmark 9.9
5 Estonia 13.1
6 Finland 12.7
7 France 6.2
8 Germany 9.7
9 Greece 8.7
10 Hungary 5.7
11 Ireland 10.4
12 Italy 8.1
13 Latvia 3.3
14 Lithuania 4.2
15 Luxembourg 24.5
16 Malta 6.3
17 Poland 8.3
18 Portugal 5.7
19 Slovakia 7.0
20 Slovenia 7.6
21 Spain 8.0
22 Sweden 5.6
23 Switzerland 5.6
24 UnitedKingdom 9.4
25 Netherlands 10.3
", sep=" ")

options(scipen=999) # suppressing scientific notation

# Name the data set and give each variable (column) a name:
emissions.data <- read.table(my.datafile, header=FALSE, col.names=c("Number", "country", "percap"))

attach(emissions.data)  # attaching the data frame

# Finding a 95% CI for the population mean per-capita emissions level, mu:

t.test(percap, conf.level = 0.95)$conf.int

# Finding a 90% CI for the population mean per-capita emissions level, mu:

t.test(percap, conf.level = 0.90)$conf.int

# Performing a t-test of H0: mu = 12 vs. Ha: mu < 12:

t.test(percap, mu=12, alternative="less")

# Performing a t-test of H0: mu = 12 vs. Ha: mu > 12:

t.test(percap, mu=12, alternative="greater")

# Performing a t-test of H0: mu = 12 vs. Ha: mu not equal 12:

t.test(percap, mu=12, alternative="two.sided")


########## If you are given only the summary statistics:

## Here, x.bar=8.82, s=4.113, n=25:

# Finding a 95% CI for the population mean per-capita emissions level, mu:

t.CI.stats(x.bar=8.82, s=4.113, n=25, conf.level=0.95)


# Performing a t-test of H0: mu = 12 vs. Ha: mu < 12:

t.test.stats(x.bar=8.82, s=4.113, n=25, mu=12, alternative="less")

# Performing a t-test of H0: mu = 12 vs. Ha: mu > 12:

t.test.stats(x.bar=8.82, s=4.113, n=25, mu=12, alternative="greater")

# Performing a t-test of H0: mu = 12 vs. Ha: mu not equal 12:

t.test.stats(x.bar=8.82, s=4.113, n=25, mu=12, alternative="two.sided")


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

## CIs for p and z-tests about p:

# Example 2:

# Reading the data into a temporary file called "my.datafile":

my.datafile <- tempfile()
cat(file=my.datafile, "  
Brocas    
Anomic    
Anomic    
Conduction
Brocas    
Conduction
Conduction
Anomic    
Conduction
Anomic    
Conduction
Brocas    
Anomic    
Brocas    
Anomic    
Anomic    
Anomic    
Conduction
Brocas    
Anomic    
Conduction
Anomic    
", sep=" ")

options(scipen=999) # suppressing scientific notation

# Name the data set and give the variable (column) a name:
aphas <- read.table(my.datafile, header=FALSE, col.names=c("type"))

# Alternatively, could type:

# aphas <- read.table("http://people.stat.sc.edu/hitchcock/aphasia.txt", header=FALSE, col.names=c("type"))

attach(aphas)  # attaching the data frame

table(type)  # Looking at a table of counts for each category

########################################################
########################################################
########################################################
#
# Copy and paste this section into R:
#
#
p.CI <- function(x, n, conf.level=0.95){
alpha <- 1-conf.level; pi.hat <- x/n
LL <- pi.hat - qnorm(1-alpha/2)*sqrt(pi.hat*(1-pi.hat)/(n))
UL <- pi.hat + qnorm(1-alpha/2)*sqrt(pi.hat*(1-pi.hat)/(n))
  my.out <- list(LL, UL, conf.level)
  names(my.out) <- c("lower limit", "upper limit", "confidence level")
return(my.out)
}
#
#
#########################################################
#########################################################
#########################################################

# Suppose p is the true proportion of aphasia sufferers having Anomic aphasia.

# Here, the number of "successes" x = 10, and the number of "trials", n = 22:

# Finding a 95% CI for the population proportion having Anomic aphasia, p:

p.CI(x=10, n=22, conf.level=0.95)

# Performing a z-test of H0: p = 0.5 vs. Ha: p < 0.5:

prop.test(x = 10, n = 22, p=0.5, alternative="less", correct=F)

# Performing a z-test of H0: p = 0.5 vs. Ha: p > 0.5:

prop.test(x = 10, n = 22, p=0.5, alternative="greater", correct=F)

# Performing a z-test of H0: p = 0.5 vs. Ha: p not equal 0.5:

prop.test(x = 10, n = 22, p=0.5, alternative="two.sided", correct=F)