#Basic Statistical Analysis in R#

#set working directory
setwd("D:/AJG/Documents/AJG Documents/R/R/RData")
getwd()

#import data
MFheights <- read.csv("MFheights.csv", header = TRUE, colClasses = c("factor", "numeric"))
str(MFheights)
summary(MFheights)
summary(MFheights$Height)

#descriptive operations
mean(MFheights$Height)
median(MFheights$Height)
sd(MFheights$Height)
var(MFheights$Height)
range(MFheights$Height) 
quantile(MFheights$Height)
IQR(MFheights$Height)
min(MFheights$Height)
max(MFheights$Height)

#simple plots
table(MFheights$Height)
boxplot(MFheights$Height)
hist(MFheights$Height)

#import data
Maths.Sport <- read.csv("maths-sport.csv", header = TRUE, colClasses = c("numeric", "numeric", "numeric"))
plot(Maths.Sport$Maths,Maths.Sport$Sport)

#correlations
cor(Maths.Sport$Maths,Maths.Sport$Sport)
cor.test(Maths.Sport$Maths,Maths.Sport$Sport)

#import data
IQ.Reading <- read.csv("IQ-Reading.csv", header = TRUE, colClasses = c("numeric", "numeric", "numeric"))
cor.test(IQ.Reading$Reading, IQ.Reading$IQ)

#regression
lm(IQ.Reading$Reading ~ IQ.Reading$IQ)

#import data
IQ.Books.Reading <- read.csv("IQ-Books-Reading.csv", header = TRUE, colClasses = c("numeric", "numeric", "numeric", "numeric"))

#multiple regression
model <- lm(IQ.Books.Reading$Reading ~ IQ.Books.Reading$IQ+ IQ.Books.Reading$Books)
summary(model)

#import data
IQ.Books.Sex.Reading <- read.csv("IQ-Books-Sex-Reading.csv", header = TRUE, colClasses = c("numeric", "numeric", "numeric", "numeric",  "numeric" ))

#multiple regression
model <- lm(IQ.Books.Sex.Reading$Reading ~ IQ.Books.Sex.Reading$IQ+ IQ.Books.Sex.Reading$Books + IQ.Books.Sex.Reading$Sex)
summary(model)

#relative importance
install.packages('QuantPsyc')
library('QuantPsyc')
lm.beta(model)

#import data
Announcement <- read.csv("Announcement.csv", header = TRUE, colClasses = c("numeric", "factor", "numeric"))
summary(Announcement)
boxplot(Announcement$Recall ~ Announcement$Method)

#t test
#t.test(numeric variable ~ dichotamous variable)
#t.test(numeric variable1, numeric variable2)
t.test(Announcement$Recall ~ Announcement$Method)

#load data
Stopping <- read.csv("alcohol-stopping.csv", header = TRUE, colClasses = c("numeric", "factor",  "numeric"))
Stopping

boxplot(Stopping$stopping ~ Stopping$alcohol)

#anova
anova1 <- aov(Stopping$stopping ~ Stopping$alcohol)
summary(anova1)
sum <-TukeyHSD(anova1)
sum
plot(sum, las=1) 
par(mar=c(5,10,5,5))

#load data
stopping <- read.csv("alcohol-gender-stopping.csv", header = TRUE, colClasses = c("numeric", "factor", "factor", "numeric"))
boxplot(stopping~alcohol, data=stopping)
boxplot(stopping~alcohol, data=stopping, col=rainbow(4))
boxplot(stopping~gender, data=stopping, col=c("red", "blue"))

#produce table of means
tapply(stopping$stopping, stopping$gender, mean)
tapply(stopping$stopping, stopping$alcohol, mean)
tapply(stopping$stopping, list(stopping$gender, stopping$alcohol), mean)

#produce table of max
tapply(stopping$stopping, stopping$gender, max)

#histogram
count <- tapply(stopping$stopping, list(stopping$gender, stopping$alcohol), mean)
barplot(count, main="Stopping Distance", beside=TRUE, xlim=c(0,16), legend=TRUE, col=c("red","blue"))

#2-way anova
anova2 = aov(stopping ~ alcohol + gender + alcohol:gender, data = stopping) 
summary(anova2)

#interaction effect
model <- lm(stopping ~ alcohol + + gender + alcohol:gender, data = stopping)
anova(model)

#interaction plot
interaction.plot(stopping$alcohol, stopping$gender, stopping$stopping)
interaction.plot(stopping$alcohol, stopping$gender, stopping$stopping, col=c("red", "blue"), trace.label="Gender", xlab="Amount of Alcohol", ylab="Stopping Distance (m)", main= "Stopping 
Distance According to Alcohol Consumption")

model_aov <- aov(model)
TukeyHSD(model_aov)

#analysis of covariance - ancova
#import data
stopping <- read.csv("alcohol-age-stopping.csv", header = TRUE, colClasses = c("numeric","factor","numeric", "numeric"))
stopping

means <- tapply(stopping$stopping, stopping$alcohol, mean)
means

ancova<-aov(stopping$stopping ~ stopping$alcohol+stopping$age)
summary(ancova)

TukeyHSD(ancova,"stopping$alcohol")

#install and load emeans package
install.packages('emmeans')
library('emmeans')

#estimate marginal means
emmeans(ancova, 'alcohol')

#two-way ancova
test <- read.csv("method-pretest-gender-test-score.csv", header = TRUE, colClasses = c("numeric", "factor","factor","numeric", "numeric"))
test
t.test(test$score ~ test$gender)

ancova2<-aov(test$score ~ test$method+test$gender+test$pretest + test$method:test$gender)
summary(ancova2)
TukeyHSD(ancova2,"test$method")
TukeyHSD(ancova2)

#install and load effects package
install.packages('effects')
library('effects')

result <- aov(score ~ pretest + method*gender, data=test)
effect_method = effect("method*gender",result)
data.frame(effect_method)

means <- data.frame(effect_method)
plot <- tapply(means$fit, list(means$gender, means$method), mean)
barplot(plot, main="Final Score, controlling for Pre-Test", beside=TRUE, xlim=c(0,20), legend=TRUE)
barplot(plot, main="Final Score, controlling for Pre-Test", beside=TRUE, xlim=c(0,20), legend=TRUE, col=c("blue","red"))
interaction.plot(means$method, means$gender, means$fit, col=c("red", "blue"), trace.label="Gender", xlab="Teaching Method", ylab="Final Score", main= "Final Score according to Teaching Method & Gender, \n controlled for Pre-Test Score")

#chi squared
smokers <- matrix(c(18,35,24,14), ncol=2)
chisq.test(smokers)

smokers <- matrix(c(18,35,24,14), ncol=2)
rownames(smokers) <- c("Male", "Female")
colnames(smokers) <- c("Smoke", "Not Smoke")
smokers
chisq.test(smokers)

#import data
smoke <- read.csv("smoke.csv", header = TRUE, colClasses = c("numeric", "factor","factor"))
smoke
table(smoke$Sex, smoke$Smoke)

MFsmoke <- table(smoke$Sex, smoke$Smoke)
MFsmoke
chisq.test(MFsmoke)