# set global options
knitr::opts_chunk$set(
  echo = TRUE, # show code 
  warning = FALSE,  
  message = FALSE
  
)

Case Study.

Productivity improvement. An economist compiled data on productivity improvements last year for a sample of firms producing electronic computing equipment. The firms were classified according to the level of their average expenditures for research and development in the past three years (low, moderate, high). The results of the study follow (productivity improvement is measured on a scale from 0 to 100 ). Assume that ANOVA model (16.2) is appropriate.

\[ \begin{array}{llccccccccccc} & & \boldsymbol{j} \\ \boldsymbol{i} & & \mathbf{1} & \mathbf{2} & \mathbf{3} & \mathbf{4} & \mathbf{5} & \mathbf{6} & \mathbf{7} & \mathbf{8} & \mathbf{9} & \mathbf{1 0} & \mathbf{1 1} & \mathbf{1 2} \\ \hline 1 & \text { Low } & 7.6 & 8.2 & 6.8 & 5.8 & 6.9 & 6.6 & 6.3 & 7.7 & 6.0 & & & \\ 2 & \text { Moderate } & 6.7 & 8.1 & 9.4 & 8.6 & 7.8 & 7.7 & 8.9 & 7.9 & 8.3 & 8.7 & 7.1 & 8.4 \\ 3 & \text { High } & 8.5 & 9.7 & 10.1 & 7.8 & 9.6 & 9.5 & & & & & & \end{array} \]

  1. Suppose a priori, the researcher wants to estimate the mean productivity improvement for firms with high research and developme expenditures levels: use a 95 percent confidence interval. Interpret your interval estimate.
Low <- c(7.6, 8.2, 6.8, 5.8, 6.9, 6.6, 6.3, 7.7, 6.0)
Moderate <- c(6.7, 8.1, 9.4, 8.6, 7.8, 7.7, 8.9, 7.9, 8.3, 8.7, 7.1, 8.4)
High <- c(8.5, 9.7, 10.1, 7.8, 9.6, 9.5)
response <- c(Low, Moderate, High)
r=3
n1=length(Low)
n2=length(Moderate)
n3=length(High)
level = c(rep(1,n1),rep(2,n2),rep(3,n3))
#means
mu_hat_overall <- mean(response)
mu_hat_low <- mean(Low)
mu_hat_moderate <- mean(Moderate)
mu_hat_high <- mean(High)
#fitted values and residuals
y_hat=c(rep(mu_hat_low,n1),rep(mu_hat_moderate,n2),rep(mu_hat_high,n3))
residuals=response-y_hat
### ANOVA Table
n_T=n1+n2+n3
SSTR = sum((y_hat-mu_hat_overall)^2)
SSE = sum((response-y_hat)^2)
SSTO = sum((response-mu_hat_overall)^2)

library(kableExtra)

AnovaTable = matrix(c(1,2,3,4,5,6,7,8,9),nrow = 3,ncol = 3)## initial setting


AnovaTable[1,1] = SSTR
AnovaTable[1,2] = r-1
AnovaTable[1,3] = AnovaTable[1,1]/AnovaTable[1,2]
 
AnovaTable[2,1] = SSE
AnovaTable[2,2] = n_T-r
AnovaTable[2,3] = AnovaTable[2,1]/AnovaTable[2,2]

SSTO = SSE + SSTR

AnovaTable[3,1] = SSTO
AnovaTable[3,2] = n_T-1
AnovaTable[3,3] = NA
 
rownames(AnovaTable) = c('Between treatments','Error (within treatments)','Total')
colnames(AnovaTable) = c('Squared Sum','Df','Mean Sum')


library(knitr)
library(kableExtra)

knitr::kable(AnovaTable, 
             "html",
             digits=6,
             caption = "ANOVA Table") %>% kable_styling()
ANOVA Table
Squared Sum Df Mean Sum
Between treatments 20.12518 2 10.062593
Error (within treatments) 15.36222 24 0.640093
Total 35.48741 26 NA 
MSE <- 0.640093

std_dev_mu_hat_high = (sqrt(.640093)/sqrt(6))


#T score obtained by T table = 2.064

CI1_upperbound <- 9.2 + 2.064*std_dev_mu_hat_high
CI1_lowerbound <- 9.2 - 2.064*std_dev_mu_hat_high

Confidence_interval_1 <- c(CI1_lowerbound, CI1_upperbound)
print("95% Confidence Interval")
## [1] "95% Confidence Interval"
Confidence_interval_1
## [1] 8.525851 9.874149

The 95% confidence interval for the factor level mean of high research and development expenditures spans 8.525851 to 9.874149. This means we may be 95% confident that the true population level factor level mean of productivity improvements for high research and development expenditures is within this interval. Or, if we were to run the experiment an infinte number of times 95% of the confidence intervals we would gather from the sample mean of the high research and development sampling distribution would contain the true population level mean.

  1. Suppose a priori, the researcher wants to estimate \(D=\mu_{2}-\mu_{1}\): use a 95 percent confidence interval. Interpret your interval estimate.
Dhat = mu_hat_moderate - mu_hat_low
std_dev_Dhat = sqrt(MSE*((1/12)+(1/9)))

#T score from T distribution table with 1-.05/2 area, and df = 24
#T score = 2.064

CI2_upperbound <- Dhat + 2.064*std_dev_Dhat
CI2_lowerbound <- Dhat - 2.064*std_dev_Dhat

Confidence_interval_2 <- c(CI2_lowerbound, CI2_upperbound)
print("95% Confidence Interval")
## [1] "95% Confidence Interval"
Confidence_interval_2
## [1] 0.5273919 1.9837192

The 95% confidence interval here spans 0.5273919 to 1.9837192. This means we are 95% confident that the true population level difference in treatment level means for the moderate and the low groups would be between 0.5273919 and 1.9837192. Alternatively, if we were to repeat this experiment infinitely many times, 95% of the time we would create a confidence interval which captures the true population level difference in treatment level means for the moderate and low groups.

  1. Suppose a priori, the researcher wants to test for all pairs of factor level means whether or not they differ. Try all appropriate procedures, which one is most efficient? Use the most efficient procedure to conduct the tests, with \(\alpha=.05 .\) Set up groups of factor levels whose means do not differ.
#Bonferroni 
g = (3*(3-1))/2
B = qt((1-(.05/(2*g))),24,lower.tail=T)

#Scheffe
S=sqrt((3-1)*qf(.95, df1=2, df2=24, lower.tail = TRUE))

#Tukey
T=(1/sqrt(2))*qtukey(.95, 3, 24, nranges = 1, lower.tail = TRUE, log.p = FALSE)

which.min(c(B,S,T))
## [1] 3
print("Use Tukey")
## [1] "Use Tukey"
#### Tukey Pair-wise Comparison

#1.) mu_hat_low - mu_hat_moderate
Dhatlowmod <- mu_hat_moderate - mu_hat_low
stddev_Dhat <- sqrt(MSE)*sqrt((1/12)+(1/9))
Test_stat = (Dhatlowmod - 0)/stddev_Dhat
if (abs(Test_stat) > T) print ("reject H0 mu_hat_low and mu_hat_moderate ARE NOT
equal")
## [1] "reject H0 mu_hat_low and mu_hat_moderate ARE NOT\nequal"
#2.)mu_hat_moderate - mu_hat_high
Dhatmodhigh <- mu_hat_high - mu_hat_moderate
stddev_Dhat <- sqrt(MSE)*sqrt((1/6)+(1/12))
Test_stat = (Dhatmodhigh - 0)/stddev_Dhat
if (abs(Test_stat) > T) print ("reject H0 mu_hat_moderate and mu_hat_high ARE 
NOT equal")
## [1] "reject H0 mu_hat_moderate and mu_hat_high ARE \nNOT equal"
#3.) mu_hat_high - mu_hat_low
Dhathighlow <- mu_hat_high - mu_hat_low
stddev_Dhat <- sqrt(MSE)*sqrt((1/6)+(1/9))
Test_stat = (Dhathighlow-0)/stddev_Dhat
if (abs(Test_stat) > T) print ("reject H0 mu_hat_high and mu_hat_low ARE 
NOT equal")
## [1] "reject H0 mu_hat_high and mu_hat_low ARE \nNOT equal"

Since we are interested in pairwise comparisons of factor level means with an established priori, we may use a Bonferroni, Sheffe, or Tukey test. To determine the most efficient we need to generate a critical value using all three tests and determine what is the smallest. Tukey generates a value of 2.497287, which is smaller than Bonferonni or Scheffe. As a result we use Tukey to perform all the individual pair-wise comparisons. These are the results Difference 1, the difference between moderate group mean minus low group mean test stastic is approximately 3.54, greater than our test statistic 2.49, reject null hypothesis that these two means are equal Difference 2, the difference between high group mean minus moderate group mean test stastic is approximately 2.68, greater than our test statistic 2.49, reject null hypothesis that these two means are equal Difference 3, the difference between high group mean minus low group mean test statistic is approximately 5.50, greater than our test statistic 2.49, reject null hypothesis that these two means are equal

  1. The sample sizes for the three factor levels are proportional to the population sizes. Suppose a priori, the economist wishes to estimate the mean productivity gain last year for all firms in the population. Estimate this overall mean productivity improvement with a 95 percent confidencc interval.
L = sum((1/3*mu_hat_low)+(1/3*mu_hat_moderate)+(1/3*mu_hat_high))
stddev_Lhat <- sqrt(MSE*(((1/3)^2)*(1/9)+((1/3)^2)*(1/12)+((1/3)^2)*(1/6)))
B=qt(1-(.05/(2*3)),df=27-3,lower.tail=TRUE)

Lowerbound = L - B*stddev_Lhat
Upperbound =  L + B*stddev_Lhat
CI = c(Lowerbound, Upperbound)
print("Confidence Interval")
## [1] "Confidence Interval"
print(CI)
## [1] 7.657923 8.482818

The economist is interested in the mean productivity gain last for all firms in the population across various factor levels. As a result we are interested in a linear combination of factor level means. Of the three family-wise multiple inference methods learned only the Bonferroni can address linear combinations.

I was confused by this question but my guess was we treat it as a linear combination where we are interested in a linear combination of each factor level mean times a coefficient of 1/3. That’s how I set up the linear the combination and the subsequent confidence interval. As such I believe the 95% confidence interval to between (7.657923,8.482818)

  1. Choose the most efficient procedure, obtain confidence intervals for the following comparisons with 90 percent family contidence coefticient: \[ \begin{array}{ll} D_{1}=\mu_{3}-\mu_{2} & D_{3}=\mu_{2}-\mu_{1} \\ D_{2}=\mu_{3}-\mu_{1} & L_{1}=\frac{\mu_{1}+\mu_{2}}{2} -\mu_{3} \end{array} \] Interpret your results and describe your findings.
S = sqrt((3-1)*qf(.90, df1=2, df2=24, lower.tail = TRUE))


#1.) Difference 1 mu3-mu2
Dhat1 <- mu_hat_high - mu_hat_moderate
Stddev_Dhat1 <- sqrt(MSE*((1/6)+(1/12)))


Lowerbound = Dhat1 - S *Stddev_Dhat1
Upperbound =  Dhat1 + S *Stddev_Dhat1
CI1 = c(Lowerbound, Upperbound)
print("Confidence Interval 1")
## [1] "Confidence Interval 1"
print(CI1)
## [1] 0.1653431 1.9679903
#2.) Difference 2 mu3-mu1
Dhat2 <- mu_hat_high - mu_hat_low
Stddev_Dhat2 <- sqrt(MSE*((1/6)+(1/9)))
Lowerbound = Dhat2 - S *Stddev_Dhat2
Upperbound =  Dhat2 + S *Stddev_Dhat2
CI2 = c(Lowerbound, Upperbound)
print("Confidence Interval 2")
## [1] "Confidence Interval 2"
print(CI2)
## [1] 1.372144 3.272301
#3 Difference 3 mu2-mu1
Dhat3 <- mu_hat_moderate - mu_hat_low
Stddev_Dhat3 <- sqrt(MSE*((1/12)+(1/9)))

Lowerbound = Dhat3 - S *Stddev_Dhat3
Upperbound =  Dhat3 + S *Stddev_Dhat3
CI3 = c(Lowerbound, Upperbound)
print("Confidence Interval 3")
## [1] "Confidence Interval 3"
print(CI3)
## [1] 0.4606629 2.0504483
#4.) Contrast between ((mu_hat_low+mu_hat_moderate)/2) - mu_hat_high
Stddev_Lhat <-  sqrt(MSE*(((.5^2)/9)+((.5^2)/12)+((1^2)/6)))
Lhat <- .5*mu_hat_low +.5*mu_hat_moderate - 1*mu_hat_high
Lowerbound  = Lhat - S*Stddev_Lhat
Upperbound = Lhat + S*Stddev_Lhat
CI4 = c(Lowerbound, Upperbound)
print("Confidence Interval 4")
## [1] "Confidence Interval 4"
print(CI4)
## [1] -2.5308374 -0.8580515

These are three differences and a contrast. Therefore as a family they are contrasts; it is not specified that the researcher has a priori. As a result the only family-wise simultaneous inference technique that can be used is Scheffe (Bonferroni requires an established priori, Tukey can only be used on pairwise differences, not contrasts) The confidence intervals are as follows: Difference 1 mu3-mu2 : (0.1653431, 1.2691970) Difference 2 mu3-mu1 : (1.372144, 3.272301) Difference 3 mu2-mu1 : (0.4606629, 2.0504483) Contrast 1 Contrast between ((mu_hat_low+mu_hat_moderate)/2) - mu_hat_high : (-2.5308374, -0.8580515)

  1. Suppose that the sample sizes have not yet been determined but it has been decided to use the same number of firms for each research and developme expenditures levels. Assume that a reasonable planning value for the error standard deviation is \(\sigma=0.5\); the range of the treatment means is \(2\); the \(\alpha\) risk is to be controlled at .01. What would be the required sample sizes if differences in the mean productivity improvement for different research and developme expenditures levels are to be detected with probability .80?
# search the smallest n value that satisfies the inequality
n=2
delta=2
sigma=.5
alpha=0.01
flag=FALSE
while(flag==FALSE){
  temp=pf(qf(1-alpha,2,3*(n-1)),df1=2,df2=3*(n-1),ncp=(n/2)*(delta/sigma)^2,lower.tail = FALSE)
  if(temp>=0.8){flag=TRUE;print(paste("The minimum sample size for each treatment:",n))}
  n=n+1
}
## [1] "The minimum sample size for each treatment: 4"

Based on the code above, the minimum sample size necessary to detect a difference between firms with the given parameters with probability of at least .80 is 4 firms per treatment .

  1. If the sample sizes determined in part (6) were employed, what would be the power of the test for treatment mean differences when \(\mu_{1}=7, \mu_{2}=8\), and \(\mu_{3}=9 ?\)
mus=c(7,8,9)
sum((mus-mean(mus))^2)*4/0.5^2
## [1] 32
power = pf(qf(1-alpha,2,3*(4-1)),df1=2,df2=3*(4-1),ncp=32,lower.tail = FALSE)

The power should be 0.9113954

  1. Suppose primary interest is in estimating the two patirwise comparisons: \[ L_{1}=\mu_{1}-\mu_{2} \quad L_{2}=\mu_{3}-\mu_{2} \] What would be the required sample sizes be if the precision of each comparison is to be \(\pm 0.5\), using the most efficient multiple comparison procedure with a 95 percent family confidence coefficient?

We don’t have a set sample size so need to try to construct a margin of error #with all three approaches to see which is the most efficent. We are assuming an established priori since this is the “primary interest” as such we can use Bonferroni. In construction of hypothetical sample sizes, Bonferroni was the most efficient. Thus, are minimum sample size per treatment is 21 samples.

print(qt(1-0.05/2,n_T-r))
## [1] 2.063899
print(sqrt((r-1)*qf(1-0.05,r-1,n_T-r)))
## [1] 2.608764
print(qtukey(1-0.05, nmeans=r, df=n_T-r)/sqrt(2))
## [1] 2.497287
#Bonferroni appears to be the most efficient

#Bonferroni

n=2
sigma=sqrt(MSE)
alpha=0.05
flag=FALSE
while(flag==FALSE){
  temp=qt(1-0.05/2,3*(n-1))*sqrt(sigma^2*(1/n+1/n))
  if(temp<=.5){flag=TRUE;print(paste("The minimum sample size for each treatment by Bonferroni:",n))}
  n=n+1
}
## [1] "The minimum sample size for each treatment by Bonferroni: 21"
#Scheffe
n=2
n_T=3*n
r=3
sigma=sqrt(MSE)
alpha=0.05
flag=FALSE
while(flag==FALSE){
  temp=(sqrt((r-1)*qf(1-0.05,r-1,n_T-r)))*sqrt(sigma^2*(1/n+1/n))
  if(temp<=.5){flag=TRUE;print(paste("The minimum sample size for each treatment by Scheffe:",n))}
  n=n+1
}
## [1] "The minimum sample size for each treatment by Scheffe: 98"
#Tukey
n=2
n_T=3*n
r=3
sigma=sqrt(MSE)
alpha=0.05
flag=FALSE
while(flag==FALSE){
  temp=(qtukey(1-0.05, nmeans=r, df=n_T-r)/sqrt(2))*sqrt(sigma^2*(1/n+1/n))
  if(temp<=.5){flag=TRUE;print(paste("The minimum sample size for each treatment by Tukey:",n))}
  n=n+1
}
## [1] "The minimum sample size for each treatment by Tukey: 90"