602-Assigment4-1

library(binom)
library(car)
library(collapsibleTree)
library(dbplyr)
library(dplyr)
library(EnvStats)
library(ggformula)
library(ggplot2)
library(gmodels)
library(htmltools)
library(ISLR)
library(knitr)
library(lawstat)
library(markdown)
library(mosaic)
library(mdsr)
library(mosaicData)
library(nycflights13)
library(olsrr)
library(plyr)
library(purrr)
library(plotly)
library(resampledata)
library(rmarkdown)
library(rpart)
library(rpart.plot)
library(rvest)
library(SDaA)
library(shiny)
library(stringi)
library(tibble)
library(tidyr)
library(tidyselect)
library(tinytex)
library(yaml)
library(shiny)

1. A pharmaceutical company conducted an experiment to compare the mean amount of time (in days) necessary for a patient to recover from the effects and complications that follow the onset of the common cold. In this experiment, thirty-one healthy females between the ages of 25 and 34 years of age were randomly divided into two groups.

In the first group, nVitC=15 females received a daily dose of 500 milligrams of Vitamin C. The second group consisted of nPlacebo=16 females, each of which was given a placebo, or rather a “fake” Vitamin C tablet. The recovery time - the time until the common cold symptoms had disappeared - for each of the thirty-one females was observed (in days). The data are provided below.

Received 500 mgs of Vitamin C: 6, 7, 7, 7, 8, 7, 7, 8, 7, 8, 10, 6, 8, 5, 6

Received a placebo: 10, 12, 8, 6, 9, 8, 11, 9, 11, 8, 12, 11, 9, 8, 10, 9

1. Do these data indicate that the recovery time is quicker with Vitamin C than without? Carry out the appropriate statistical investigation with a permutation test. Ensure you show “where” your permutation test statistic lies on the distribution. Report your (empirical) P-value and your statistical inference. For the number of permutation tests, use 2999.

Prepare data, display boxplot, run basic stats

time = c(6, 7, 7, 7, 8, 7, 7, 8, 7, 8, 10, 6, 8, 5, 6, 10, 12, 8, 6, 9, 8, 11, 9, 11, 8, 12, 11, 9, 8, 10, 9)
ttment = c('vc','vc','vc','vc','vc','vc','vc','vc','vc','vc','vc','vc','vc','vc','vc', 'n','n','n','n','n','n','n','n','n','n','n','n','n','n','n','n')

df = data.frame(time, ttment)
head(df, 3)

ggplot(df, aes(x=ttment, y=time)) + geom_boxplot(col='red', fill='blue', width=0.5) + ggtitle('Cold recovery time comparing treatment vs no treatment') + coord_flip()

stats = favstats(~time | ttment, data = df)
stats

State Hypotheses H0: recovery time without treatment = recovery time with treatment –> recovery time without treatment - recovery time with treatment = 0 HA: recovery time without treatment != recovery time with treatment –> recovery time without treatment - recovery time with treatment > 0

Possible permutations of 31 people choosing 15 at the time

choose(31, 15)

[1] 300540195

proceed to do a random distribution of difference of means by conducting a permutation test by sampling 15 data point without replacement as treatment and the other 16 to the no treatment, calculate the means difference between the groups and store this difference into a vector for further calculations after we replicate this process 4999 more times as the first difference in means correspond to the observed values. Total permutations 5000

Determine the mean difference between the groups

ttmentMeanDiff = stats$mean[1] - stats$mean[2]
ttmentMeanDiff

[1] 2.304167

N = 2999
meanDiff.Per = numeric(N) #create a vector to store differences of means
index = sample(31, 15, replace=FALSE) #randomly pick 15 out of 31
index

 [1] 13 18  9  5 19 24 15 31 29 20 27 28 26 16  7

for(i in 1:N)
{ index = sample(31, 15, replace=FALSE) #randomly pick 15 numbers out of 31
  meanDiff.Per[i] = mean(df$time[index]) - mean(df$time[-index]) #mean difference between recovery time no treatment and treatment
}

Now we convert the meanDiff vecto into a dataframe

perm.df = data.frame(meanDiff.Per)
head(perm.df, 4)

tail(perm.df, 4)

Display the distribution of man difference and the value for the observed difference

hist(meanDiff.Per, xlab='Difference in recovery time between treated and untreated patients', ylab='Frequency', main='Permutation Test (3000) outcome', col='darkgreen')
abline(v=ttmentMeanDiff, col='red')

Compute Empirical P-value

pval = (1 + sum(meanDiff.Per >= (ttmentMeanDiff))) / (N + 1)
pval

[1] 0.0003333333

Getting a Pvalue of 0.00033 leads us to reject the hypothesis null indicating that the observed difference in means is larger than expected if it was just random plus the observed difference in recovery times is greater the values obtained in the permutation test. In conclusion we rehect H0 and conclude that there is a relation between treatment with Vit C and faster recovery times.

t.test(~time|ttment, alternative="greater", data=df)

    Welch Two Sample t-test

data:  time by ttment
t = 4.445, df = 27.08, p-value = 6.723e-05
alternative hypothesis: true difference in means between group n and group vc is greater than 0
95 percent confidence interval:
 1.421325      Inf
sample estimates:
 mean in group n mean in group vc 
        9.437500         7.133333 

Pvalue from T test is smaller than Empirical Pvalue.

2. Re-test your statistical hypothesis in part (a) using the t-test. In doing so, state any assumptions about these data or conditions you are imposing on these data and conduct the necessary diagnostics to either confirm or refute such assumptions. Ensure you provide both the P-value and its interpretation related to these data.

We assumed the observed data (both Vit C treatment and placebo) have a normal distribution. We will display the Normal probability plot to confirm this assumption on both Vitamin c treatment and Placebo groups of observed values

ggplot(filter(df, ttment=="vc"), aes(sample= time)) + stat_qq(col='blue') + stat_qqline(col='red') + ggtitle('Normal Probability Distribution of Cold recovery Time: Vitamin c treatment')

ggplot(filter(df, ttment=="n"), aes(sample= time)) + stat_qq(col='blue') + stat_qqline(col='red') + ggtitle('Normal Probability Distribution of Cold recovery Time: Placebo')

2. In a 1981 case , McCleskey v. Zant, lawyers for the defendent demonstrated that amongst all African-American convicted murderers in the state of Georgia, 35% were found guilty of killing caucasians and 6% of those who killed African-Americans received the death penalty. As a result, the lawyers claimed that racial discrimination was a factor in sentencing. This was contested by lawyers for the State of Georgia. They argued that “murders of black victims were more likely to be unaggravated bar-room brawls, liquor-induced arguments, or lovers quarrels” - crimes that rarely result in a death sentence. But, amongst white-victims, there was a higher proportion of killings committed in the course of an armed robbery or involving torture - which more often result in death sentences.

The table below summarizes data for all African-American convicted murders in Georgia, resulting from difference categories of crime. (For example, Aggravation Level 1 consists of bar-room brawls, liquor-induced arguments, etc., Aggravation Level 6 includes the most vicious, cold-blooded, unprovoked crimes).

                                          Sentenced to Death        Not Sentenced to Death
                 v  Victim was Caucasian          45                          85
             Victim was African-American          14                         218
             

1. Visualize these data with a bar-graph (see Assignment One) that appropriately summarizes these counts.

conviction = rbind(c(45, 85), c(14, 218))
conviction

     [,1] [,2]
[1,]   45   85
[2,]   14  218

rownames(conviction) = c("Victim was Caucasian", "Victim as African American")
colnames(conviction) = c("Sentenced to Death", "Not Sentenced to Death")
conviction

                           Sentenced to Death Not Sentenced to Death
Victim was Caucasian                       45                     85
Victim as African American                 14                    218

barplot(conviction, main="Convictions based on the  victim's race", col=c("darkblue", "darkred"), legend=rownames(conviction)) 

2. Do these data suggest that the race of the victim does appear to affect whether an African-American convicted of murder in Georgia will receive a death sentence? Interpret your findings. Should you reject the null hypothesis, provide an interval with 95% coverage/confidence.

Hypotheses H0: The race of the victim and the sentence are independent HA: the race of the victim and the sentence are not independent

xchisq.test(conviction, correct=FALSE)

    Pearson's Chi-squared test

data:  x
X-squared = 49.888, df = 1, p-value = 1.628e-12

  45.00    85.00 
( 21.19) (108.81)
[26.76]  [ 5.21] 
< 5.17>  <-2.28> 
   
  14.00   218.00 
( 37.81) (194.19)
[15.00]  [ 2.92] 
<-3.87>  < 1.71> 
   
key:
    observed
    (expected)
    [contribution to X-squared]
    <Pearson residual>

prop.test(c(45, 14), c(130, 232), alternative="greater", correct=FALSE)

    2-sample test for equality of proportions without continuity correction

data:  c out of c45 out of 13014 out of 232
X-squared = 49.888, df = 1, p-value = 8.138e-13
alternative hypothesis: greater
95 percent confidence interval:
 0.2125175 1.0000000
sample estimates:
    prop 1     prop 2 
0.34615385 0.06034483 

through chisq test, we get a Pvalue that is very small, 8.138e-13, which is confirmed through the prop test with Pvalue of 8.14e-13; the confidence interval found through the prop test gives us 0.2125 - 1. We got a confidence interval betweeen 0.2125 and 1. These results leads us to reject the H) and conclude that there is a relation between the race of the victim and the sentence.

Zobs = sqrt(49.888)
Zobs

[1] 7.063144

3. Refer to a question from Assignment 1, where you were asked to refer to certain variables of the General Society Survey of 2002. For your convenience, the data file is linked below.

gss = read.csv("https://raw.githubusercontent.com/Statman44/Data602/main/GSS2002.csv")
head(gss, 3)

gss.df = gss[complete.cases(gss),]
gss.df

1. Is there a relationship between one’s support for gun laws (variable name is GunLaw) and their opinion about current government spending on Science (variable name is SpendSci)? State the appropriate statistical hypotheses.

tableGss = tally(~GunLaw + SpendSci, data=gss.df)
tableGss

        SpendSci
GunLaw   About right Too little Too much
  Favor           39         24        6
  Oppose           9         11        3

Hypotheses H0: Gun Law and Spending on Science support are Independent of each other HA: Gun Law and Spending on Science support are not independent of each other

2. Carry out the appropriate statistical test, providing both the test statistic and the P-value.

xchisq.test(tableGss, simulate.p.value=TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  x
X-squared = 2.1048, df = NA, p-value = 0.3793

   39       24        6   
(36.00)  (26.25)  ( 6.75) 
[0.250]  [0.193]  [0.083] 
< 0.50>  <-0.44>  <-0.29> 
     
    9       11        3   
(12.00)  ( 8.75)  ( 2.25) 
[0.750]  [0.579]  [0.250] 
<-0.87>  < 0.76>  < 0.50> 
     
key:
    observed
    (expected)
    [contribution to X-squared]
    <Pearson residual>

According to the chi squared test we get a Pvalue of 0.918 and Chi-squared of 2.1048

3. What can you conclude? Do these data support your null hypothesis in part (a)? State your decision and conclusion. According to the results we got, we fail to reject the hypothesis null. We can say that there is no relation between the support of Gun Law and Spending on Science.

4. Re-trace your result, in the form of a bar-graph, that was provided in Assignment 1, Question 9. Can you infer from these data that one’s level of Education is independent of their Race? Present your findings in the form of a paragraph, outlining the decision you have made, why you made the decision you made, and the P -value.

ggplot(data=gss) + geom_bar(aes(x = Education, fill = Race), position = "fill")

ggplot(data=gss) + geom_bar(aes(x = Education, fill=Race), position = "dodge", na.rm=TRUE)

tableGss1 = tally(~Education + Race, data=gss.df)
tableGss1

           Race
Education   Black Other White
  Bachelors     2     1    24
  Graduate      1     1    15
  HS            5     2    29
  Jr Col        1     0     7
  Left HS       0     0     4

Hypotheses: H0: the level of education and Race are independent HA: the level of education and race are not independent

xchisq.test(tableGss1, simulate.p.value=TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  x
X-squared = 2.5062, df = NA, p-value = 0.979

    2        1       24   
( 2.64)  ( 1.17)  (23.18) 
[0.1557] [0.0258] [0.0287]
<-0.39>  <-0.16>  < 0.17> 
     
    1        1       15   
( 1.66)  ( 0.74)  (14.60) 
[0.2644] [0.0921] [0.0111]
<-0.51>  < 0.30>  < 0.11> 
     
    5        2       29   
( 3.52)  ( 1.57)  (30.91) 
[0.6205] [0.1208] [0.1184]
< 0.79>  < 0.35>  <-0.34> 
     
    1        0        7   
( 0.78)  ( 0.35)  ( 6.87) 
[0.0604] [0.3478] [0.0025]
< 0.25>  <-0.59>  < 0.05> 
     
    0        0        4   
( 0.39)  ( 0.17)  ( 3.43) 
[0.3913] [0.1739] [0.0930]
<-0.63>  <-0.42>  < 0.30> 
     
key:
    observed
    (expected)
    [contribution to X-squared]
    <Pearson residual>

According to the results, Pvalue=0.759 so we fail to reject the hypothesis null and conclude that there is no relation between one’s level of education and race.

4. A group of patients with a binge-eating disorder were randomly assigned to take either the experimental drug fluvoxamine or the placebo in a nine-week-long, double-blinded clinical trial. At the end of the trial the condition of each patient was classified into one of four categories: no response, moderate response, marked response, or remission. The table below shows a cross-classification, or contingency table, of these data.

          No Response        Moderate Response       Marked response     Remission     

   Fluvoxamine 15 7 3 15 Placebo 22 7 3 11

Do these data provide statistically significant evidence to conclude that there is an association between the type of treatment received and a patient’s response?

Hypotheses H0: The patient’s response to treatment and the use of Fluvoxamine are independent HA: The patient’s response to treatment and the use of Fluvoxamine are not independent

bingeTable = rbind(c(15, 7, 3, 15), c(22, 7, 3, 11))
rownames(bingeTable) = c("Fluvoxamine", "Placebo")
colnames(bingeTable) = c("No Response", "Moderate Response", "Marked Response", "Remission")
bingeTable

            No Response Moderate Response Marked Response Remission
Fluvoxamine          15                 7               3        15
Placebo              22                 7               3        11

xchisq.test(bingeTable, simulate.p.value=TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  x
X-squared = 1.8337, df = NA, p-value = 0.6332

 15.00     7.00     3.00    15.00  
(17.83)  ( 6.75)  ( 2.89)  (12.53) 
[0.4496] [0.0095] [0.0041] [0.4869]
<-0.670> < 0.097> < 0.064> < 0.698>
       
 22.00     7.00     3.00    11.00  
(19.17)  ( 7.25)  ( 3.11)  (13.47) 
[0.4182] [0.0088] [0.0038] [0.4529]
< 0.647> <-0.094> <-0.062> <-0.673>
       
key:
    observed
    (expected)
    [contribution to X-squared]
    <Pearson residual>

We got a Pvalue of 0.6082 we fail to reject the hypothesis null and therefore we can say that there is no relation between the treatment with Fluvoxamine and the response of the patients.

5. A study conducted by Youmans and Jee2 looked at students who took a research methods class in psychology. Two lecture sections of the same course were taught by the same instructor. Students in each lecture section were required to register in one of two Friday afternoon discussion sections, or a tutorial. During the ninth week of instruction, all students were asked to fill out an informal midsemester evaluation in their Friday afternoon discussion session.

In half of the Friday afternoon discussion sections, the experimenter (not the professor/instructor!) passed around a bag containing small bars of chocolate that he simply had ‘left-over’ and ‘wanted to get rid of’ so that students would not think the chocolate was a gift from their professor/instructor. The evaluation was a survey that contained nine questions. For each question the student provided a rating from 1 (very poor) to 5 (excellent).

Question 9 posed the following: ``Compared to other instructors you have had at this university, this instrutor is:’’ 1 (very poor) to 5 (excellent)

Compliments to the authors, they have provided me with the data in this particular study. In the file you will find the raw data providing the student responses for Question 9 (Q9) and the overall average response for Questions 1 through 9 (Overall) for each student. Note: To easily read this into your R session, select the “Raw” button prior, then use the resulting URL.

Ssurvey = read.csv("https://raw.githubusercontent.com/Statman44/Data602/main/chocnochocratings.csv")
head(Ssurvey, 3)

1. Pertaining to Question 9: Do these data suggest there is a treatment effect? Test the existence of a treatment effect using a permutation test. If a treatment effect is discovered, explain its meaning in the context of these data. Carry out a permutation test, using 1999 iterations.

(You will have to carefully consider “what a treatment effect” is in this scenario.)

The treatment here is the administration of chocolate. By looking at the data we can see that they all were offered the chocolate but some took it some others didn’t take the chocolate. Did the administration of chocolate influence their survey response to Q9?

Hypotheses H0: The responses to Q9 in the survey and the treatment are independent Q9 treatment = Q9 no treatment HA: the responses to Q9 in the survey and the treatment are not independent Q9 treatment > Q9 no treatment

Exploring the data

ggplot(Ssurvey, aes(x = GroupName, y = Q9)) + geom_boxplot(col="red", fill="blue", width=0.5) + ggtitle("Q9 Student Survey: treatment effect") + coord_flip()

favstats(~Q9|GroupName, data= Ssurvey)

obsdiff.ttmt = mean(~Q9, data = filter(Ssurvey, GroupName == "Chocolate")) - mean(~Q9, data = filter(Ssurvey, GroupName == "NOChoc"))
obsdiff.ttmt

[1] 0.2033333

Performing 1999 permutations

N = 1999 
perms = numeric(N) 
for(i in 1:N)
{ index = sample(98, 50, replace=FALSE) 
  perms[i] = mean(Ssurvey$Q9[index]) - mean(Ssurvey$Q9[-index]) 
}

perms.df = data.frame(perms)
head(perms.df,4)

tail(perms.df, 3)

Display a histogram to understand the distribution of the permutations

hist(perms, xlab="Difference Between Mean answer to Q9 (Chocolate - NoChocolate)", ylab="Frequency", main="Outcome of 2000 Permutation Tests", col='blue')
abline(v = obsdiff.ttmt, col="red")

Find the Empirical Pvalue

(sum(perms >= obsdiff.ttmt) + 1 )/(N + 1)  #computes empirical P-value

[1] 0.161

2. Consider the variable Overall. Is there a treatment effect with respect to the professor’s overall rating as a teacher? Apply the t-test to these data. Interpret the meaning of the P-value. (chitest and xchitest)

Hypotheses H0: the overall score and the treatment are independent HA: the overall score and the treatment are not independent

tableSsurvey = tally(~Overall + GroupName, data=Ssurvey)
tableSsurvey

       GroupName
Overall Chocolate NOChoc
   2            0      1
   2.44         1      0
   2.56         1      1
   2.67         1      0
   2.78         0      1
   2.89         2      1
   3            2      3
   3.11         0      2
   3.33         0      2
   3.44         1      1
   3.56         1      2
   3.67         6      2
   3.78         2      4
   3.89         2      6
   4            3      3
   4.11         3      0
   4.22         3      8
   4.33         2      2
   4.44         3      2
   4.56         6      3
   4.67         3      2
   4.78         1      0
   4.89         3      0
   5            4      2

Graphing Normal Probability plot for both Treatment and No Treatment

ggplot(filter(Ssurvey, GroupName == "Chocolate"), aes(sample = Overall)) + stat_qq(col="blue") + stat_qqline(col="red") + ggtitle("Normal Probability Plot of Student Survey: Treatment")

ggplot(filter(Ssurvey, GroupName == "NOChoc"), aes(sample = Overall)) + stat_qq(col="blue") + stat_qqline(col="red") + ggtitle("Normal Probability Plot of Student Survey: No Treatment")

t.test(~Overall | GroupName, data = Ssurvey, alternative='greater')

    Welch Two Sample t-test

data:  Overall by GroupName
t = 1.6616, df = 95.93, p-value = 0.04993
alternative hypothesis: true difference in means between group Chocolate and group NOChoc is greater than 0
95 percent confidence interval:
 9.400666e-05          Inf
sample estimates:
mean in group Chocolate    mean in group NOChoc 
               4.072000                3.849792 

According the the Pvalue we get from the t.test of 0.05, and the empirical Pvalue of 0.161, we fail to reject the hypothesis null and conclude that the professor’s overall rating is independent of the student’s treatment with chocolate

chisq.test(tableSsurvey, simulate.p.value=TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  tableSsurvey
X-squared = 24.842, df = NA, p-value = 0.3708

3. Consider the test suggested in part (a). Why would a t-test not be a recommended statistical method to carry out the test in part (a)? Explain your answer in a few sentences.

R./ We can see in the boxplots how similar the stats for the 2 groups are, this is confirmed by the favstats results. The permutation will help magnify the differences between the two groups and give statistical significance to the difference, if any. Permutations are good when we do not know much about the population and we don’t need to do assumptions. The Ttest assumes the data has normal distribution whiche here we have some data skewed.