Enter your search text in the box above

Select a result to preview

DiscoveringStatsUsingR

本书章节对应年级（英制）

Chapter 1st year undergraduate 2nd year undergraduate Post-graduates
1 ✓
2 ✓
3 ✓
4 ✓
5 ✓ ✓
6 ✓
7 ✓
8 ✓
9 ✓
10 ✓
11 ✓
12 ✓
13 ✓
14 ✓
15 ✓
16 ✓
17 ✓
18 ✓
19 ✓
Levels ① ② ③ ③ ④

BODMAS 运算顺序

C1. Why to learn statistics?

从提出研究问题到分析数据

graph TD

A[Research questions]--> B1[data]
A[Research questions]--> B2[提出theory]

B2 -->|需要quantify| C[提出hypothesis]
C -->|需要确定测量的变量vars| C1[收数据、检验theory]
C --> C2[identify vars]
	C1 --> D[分析数据]
	C1 --> E[measure vars]
	C2 --> C3(自变量，因变量)
D --> B2
D ==> F[graph data; fit a model]



Z[data贯穿整个研究]

科学理论需要以下属性

1.5.1.1. Indepen. Var & depen. Var

1.5.1.2. Levels of measurement

Categorical var
- Human
A bit of human (×)
- Male/female → Binary var (只可分为 2 类)
- 冠/亚/季军
这些变量间不可进行数学运算；如：冠军 × 亚军=啥？！		 → nominal var (分类≥3)
- 冠/亚/季军
变量间有次序之分，但还是不能量化变量间的差异		 → 同时也是 ordinal var
Continuous var
- 量表
- 5-point scale; 不同 point 之间的 interval 是相等的		 → Interval var

Zero is arbitrary
0 也意味着东西，如温度 0
- RT
100ms VS 200ms
200ms VS 400ms
200ms 是 100ms 的 2 倍，是 400ms 的 1/2		 → Ratio var

Zero is absolute
0 就是无，如时间
- 可分为 continuous 和 discrete
- Continuous 如年龄可以是 34 yrs 7 mths 21 days 55mins 10 s
- Discrete 常为整数

1.5.3. Validity & reliability

1.6.2. Experimental research methods

1.6.2.1. Two methods of data collection

1.6.2.2. Two types of variation

1.7. Analyzing data

1.7.1. Frequency distributions

Pasted image 20240709105659.png

1.7.2. The center of a distribution

1.7.2.1. The mode (P22)

1.7.2.2. The median

1.7.2.3. The mean

1.7.3. The dispersion in a distribution

1.7.4. Using a frequency distribution to go beyond the data

z-score probability
±1.96 之间 95%
±2.58 之间 99%
±3.29 之间 99.9%

Pasted image 20240928162432.png

1.7.5. Fitting statistical model to the data

SD VS SE

总结来说，标准差衡量数据本身的变异性，而标准误衡量样本均值对总体均值的估计精度。随着样本量增加，标准误会减小，表示估计的准确性提高。

C2. Everything about statistics

2.2. Building statistical models (P34)

2.4. Simple statistical models (P36)

2.4.1. The mean: a very simple stats model

2.4.2. Assessing the fit of the mean: sums of squares 平方和, variance 方差, SD 标准差

如何理解自由度 df？

2.4.3. Expressing the mean as a model (P40)

2.5. Going beyond the data

2.5.1. The standard error (SE) 标准误

2.5.2. Confidence intervals 置信区间

2.5.2.1. Calculating CI

2.5.2.2. Calculating other confidence intervals

For small sample size

We’re interested to know how many followers people have with on Instagram. We collect a sample of 11 individuals and then calculate the mean and sd. We find that the mean is 96.64 and the standard deviation is 61.27.
a. Calculate a 95% confidence interval for this mean.
Since the sample size is 11, less than 30, we need to use the t-score to calculate the CI

# input data
follower_mean <- 96.64       # Sample mean
follower_sd <- 61.27         # Sample standard deviation
n <- 11                 # Sample size
confidence_level <- 0.95  # 95% confidence level

# Step 1: Calculate the degrees of freedom
df <- n - 1

# Step 2: Get the critical t-score
t_score <- qt(1 - (1 - confidence_level) / 2, df) # that returns the quantile (or critical value) from the t-distribution for a specified probability and degrees of freedom.

# Step 3: Calculate the standard error
standard_error <- follower_sd / sqrt(n)

# Step 4: Calculate the margin of error
margin_of_error <- t_score * standard_error

# Step 5: Calculate the confidence interval
lower_bound <- follower_mean - margin_of_error
upper_bound <- follower_mean + margin_of_error

# Display the results
cat("Sample Mean:", follower_mean, "\n")
cat("T-score for", confidence_level * 100, "% confidence level with", df, "df:", round(t_score, 3), "\n")
cat("Margin of Error:", round(margin_of_error, 3), "\n")
cat("Confidence Interval: [", round(lower_bound, 3), ", ", round(upper_bound, 3), "]\n")
For large sample size

b. We collect new data and now we have 56 subject total. Recalculate the confidence interval.
Now the sample size is over 30, we need to use z-score to calculate the CI

# Step 1: Calculate sample statistics
follower2_mean <- 96.64  # Sample mean
follower2_mean
follower2_sd <- 61.27                      # Population standard deviation (known)
follower2_sd
n <- 56          # Sample size
n

# Step 2: Set confidence level and find the z-score
confidence_level <- 0.95          # 95% confidence level
alpha <- 1 - confidence_level     # Significance level (0.05)
z_score <- qnorm(1 - alpha / 2)   # Z-score for 95% confidence. qnorm() returns the quantile (or critical value) of the normal distribution for a given cumulative probability. its like feeding it qnorm(0.975)
z_score

# Step 3: Calculate the standard error
standard_error <- follower2_sd / sqrt(n)  # Standard Error

# Step 4: Calculate the margin of error and confidence interval
margin_of_error <- z_score * standard_error 
lower_bound <- follower2_mean - margin_of_error
upper_bound <- follower2_mean + margin_of_error

# Display the results
cat("Sample Mean:", round(follower2_mean, 2), "\n")
cat("Z-score for", confidence_level * 100, "% confidence level:", round(z_score, 2), "\n")
cat("Margin of Error:", round(margin_of_error, 2), "\n")
cat("Confidence Interval: [", round(lower_bound, 2), ",", round(upper_bound, 2), "]\n")

2.5.2.4. Showing confidence intervals visually

2.6. Using statistical models to test research questions

2.6.1. Test statistics

P(1∪2∪3)=16+16+16=12

Jane superbrain 2.6 (P54)

CAN CANNOT
1 "there is a significant effect of ..." "the effect is important"

因为，very small and unimportant effects can turn out to be statistically significant just because huge numbers of people have been used
只要招的人够多，再小、再不重要的效应都可以变成统计学意义上的“显著”效应
2 若算出来的 test statistic occurring by chance 的概率 > 0.05, 我们可以 reject H1 "the H0 is true"

因为 H0 的含义是：there is no effect in the population；就算是 test statistic > 0.05, 也只能说明 H0 is true 的概率很大，而不是说 H0 就是true
3 若算出来的 test statistic occurring by chance 的概率 < 0.05, 我们可以 reject H0 "the H0 is false"

同理，只能说明 H0 is false 的概率更大，但不能说明 H0 就是false

2.6.2. One-and two-tailed tests

url: https://www.bilibili.com/video/BV1oP4y177Qh/?spm_id_from=333.788&vd_source=f1111db8c0296ebbc68186619f534b4f
title: "通俗统计学原理入门7 单边检验 单边假设检验 单边t检验 单尾检验 单侧检验_哔哩哔哩_bilibili"
description: "新年，收获了人生的第一个1000个粉丝。觉得很好玩，完全是无心插柳。本来是带着功利心开始做这个统计学微课视频的，没想到无意中帮助了很多陌生人。感觉冥冥中，印证了“应无所住，而生其心”这句话。当然，这句话也是我对正态分布和中心极限定理的领悟。本学期我也确实上了这门课，见证过同学们“领悟”的那一刻。他们“哦”的那个表情，是不会骗人的。那个瞬间，我是最欣慰的，感觉像度了一个人一样。粉丝留言里，还真有一位, 视频播放量 84689、弹幕量 340、点赞数 2241、投硬币枚数 1722、收藏人数 1773、转发人数 506, 视频作者 陈祥雨大猫咪老师, 作者简介 ，相关视频：SPSS非参数检验视频教程汇总（非常重要，建议收藏！），通俗统计学原理入门12 - 置信区间 Confidence Interval 区间估计 点估计 t临界值 标准误，通俗统计学原理入门10 p值的含义 p value，通俗统计学原理入门6 关键一步 从均值抽样分布到t分布，一招搞定单/双侧检验判断｜假设检验｜统计学｜快速口诀，通俗统计学原理入门8 自由度 单样本t检验自由度，从图上看假设检验里的两类错误以及它们的关系，如何选择统计学方法？T检验、单因素方差分析、秩和检验、卡方检验到底应该选择哪一个？一个视频轻松搞定，通俗统计学原理入门 - t检验 正态分布 显著水平，关于假设检验的一切 - 统计学"
host: www.bilibili.com
image: //i0.hdslb.com/bfs/archive/71af0de10c4d4d71c65398c560994850f0d2f677.jpg@100w_100h_1c.png

Pasted image 20241004120309.png

Pasted image 20241004120336.png

Pasted image 20241004120437.png

Pasted image 20241004120506.png

Pasted image 20241004120539.png

Pasted image 20241004120605.png

Pasted image 20241004120629.png

∵原假设从双边变成了单边
Pasted image 20241004120703.png

2.6.3. Type I & Type II errors

实际上：被告无罪（Null 为真） 实际上被告有罪（Null 为假）
判无罪 ✅ 正确 ❌ Type II Error（漏判）
判有罪 ❌ Type I Error（误判） ✅ 正确

Pasted image 20241108155319.png

Pasted image 20241108155246.png

2.6.4. Effect size

Meta-analysis

2.6.5. Statistical power

C3. The R environment

见 E:/R/dsur/C3R_environment/Chapter 3 The R Environment.R

C4. Exploring Data with Graphs

url: https://www.color-hex.com/
title: "Color Hex Color Codes"
description: "Color hex is a easy to use tool to get the color codes information including color models (RGB,HSL,HSV and CMYK), css and html color codes."
host: www.color-hex.com
favicon: /favicon.ico

4.2.2. What makes a good graph (7 points)

画图应注意的点

4.4. Introducing ggplot2

4.4.1. The anatomy of a plot

Optional_页面_2.jpg

4.4.3. Aesthetics

4.4.5. Stats and geoms

4.4.6. Avoiding overplotting

4.4.7. Saving graphs

ggsave(filename)
# e.g.
ggsave("outlier amazon.png", width = 2, height = 2)

4.4.8. Quick tutorial

#set cran for r code running in obsidian 
options(repos = c(CRAN = "https://cran.r-project.org"))

#Set the working directory (you will need to edit this to be the directory where you have stored the data files for this Chapter)

setwd("E:/R/dsur/C4exploring_data_with_graphs")
imageDirectory<-"E:/R/dsur/C4exploring_data_with_graphs/C4_images"

#imageDirectory<-file.path(Sys.getenv("HOME"), "Documents", "Academic", "Books", "Discovering Statistics", "DSU R", "DSU R I", "DSUR I Images")


######A function to make it quick to save graphs in the image directory

saveInImageDirectory<-function(filename){
	imageFile <- file.path(imageDirectory, filename)
	ggsave(imageFile)	
}



######Initiate packages

#If you don't have ggplot2 installed then use:
install.packages(c("ggplot2", "plyr"))


#Initiate ggplot2
library(ggplot2)
library(reshape)
library(plyr)


#--------4.4.8.Quick Tutorial----------
facebookData <- read.delim("FacebookNarcissism.dat.txt",  header = TRUE)

graph <- ggplot(facebookData, aes(x = NPQC_R_Total, y = Rating))
graph + 
    geom_point() + 
    ggtitle("geom_point()")

Pasted image 20241031172742.png

graph + 
    geom_point(shape = 17) + 
    ggtitle("geom_point(shape = 17)") #shape = 17 is a triangle

Pasted image 20241031172757.png

graph + 
    geom_point(size = 6, shape = 17) + 
    ggtitle("geom_point(size = 6)")

Pasted image 20241031172812.png

graph + 
    geom_point(aes(colour = Rating_Type)) + 
    ggtitle("geom_point(aes(colour = Rating_Type))")

Pasted image 20241031172842.png

graph + 
    geom_point(aes(colour = Rating_Type), position = "jitter") + 
    ggtitle("geom_point(aes(colour = Rating_Type), position = jitter)")

Pasted image 20241031172858.png

graph + 
    geom_point(aes(shape = Rating_Type), position = "jitter") + 
    ggtitle("geom_point(aes(shape = Rating_Type), position = jitter)")

Pasted image 20241031173112.png

4.5. Scatterplot

4.5.1. Simple scatterplot

examData <- read.delim("Exam Anxiety.dat.txt",  header = TRUE)
head(examData)
names(examData)

# Simple scatter
scatter <- ggplot(examData, aes(x = Anxiety, y = Exam))
scatter + 
    geom_point() + 
    labs(x = "Exam Anxiety", y = "Exam Performance %")

Pasted image 20241031173134.png

4.5.2. Adding a funky line

#Simple scatter with geom_smooth()
scatter <- ggplot(examData, aes(x = Anxiety, y = Exam))
scatter + 
    geom_point() + 
    geom_smooth() + 
    labs(x = "Exam Anxiety", y = "Exam Performance %")

Pasted image 20241031173201.png

#Simple scatter with regression line
scatter <- ggplot(examData, aes(x = Anxiety, y = Exam))
scatter + 
    geom_point() + 
    geom_smooth(method = "lm", colour = "Red") + 
    labs(x = "Exam Anxiety", y = "Exam Performance %")

Pasted image 20241031174001.png

## "se = F" means "standard error = False", can switch off the confidence interval
scatter <- ggplot(examData, aes(x = Anxiety, y = Exam))
scatter + 
    geom_point() + 
    geom_smooth(method = "lm", colour = "Red", se = F) + 
    labs(x = "Exam Anxiety", y = "Exam Performance %")

Pasted image 20241031174040.png

#Simple scatter with regression line + colored CI
scatter <- ggplot(examData, aes(x = Anxiety, y = Exam))
scatter + 
    geom_point() + 
    geom_smooth(method = "lm", colour = "Red", alpha = 0.1, fill = "Red") + 
    labs(x = "Exam Anxiety", y = "Exam Performance %")

Pasted image 20241031174217.png

4.5.3. Grouped scatterplot

#Simple grouped scatter with regression line
scatter <- ggplot(examData, aes(x = Anxiety, y = Exam, color = Gender))
scatter + 
    geom_point() + 
    geom_smooth(method = "lm", colour = "Red") + 
    labs(x = "Exam Anxiety", y = "Exam Performance %")

Pasted image 20241031174329.png

scatter <- ggplot(examData, aes(x = Anxiety, y = Exam, color = Gender))
scatter + 
    geom_point() + 
    geom_smooth(method = "lm", aes(fill = Gender), alpha = 0.1)

Pasted image 20241031174650.png

scatter <- ggplot(examData, aes(x = Anxiety, y = Exam, color = Gender))
scatter + 
	geom_point() + 
	geom_smooth(method = "lm", aes(fill = Gender), alpha = 0.1) +
	labs(x = "Exam Anxiety", y = "Exam Performance %", color = "Gender")

Pasted image 20241031174845.png
（跟上面的一样诶）

scatter <- ggplot(examData, aes(x = Anxiety, y = Exam, color = Gender))
scatter + 
    geom_point() + 
    geom_smooth(method = "lm", aes(fill = Gender), alpha = 0.1) + 
    scale_fill_manual(values = c("#22ff22", "yellow")) +
    labs(x = "Exam Anxiety", y = "Exam Performance %", color = "Gender")

Pasted image 20241031174918.png

4.6. Histograms: spot obvious problems

#导入数据
festivalData <- read.delim("DownloadFestival.dat.txt",  header = TRUE)

festivalHistogram <- ggplot(festivalData, aes(day1)) + theme(legend.position="none")
festivalHistogram + 
    geom_histogram(binwidth = 0.4) + 
    labs(x = "Hygiene (Day 1 of Festival)", y = "Frequency")

Pasted image 20241028214225.png

festivalData2 = read.delim("DownloadFestival(No Outlier).dat.txt",  header = TRUE)

festivalDensity <- ggplot(festivalData2, aes(day1))
festivalDensity + 
    geom_density(aes(fill = gender), alpha = 0.5) + 
    labs(x = "Hygiene (Day 1 of Festival)", y = "Density Estimate")

Pasted image 20241028215025.png

festivalDensity + 
    geom_histogram(binwidth = 0.4, aes(fill = gender)) + 
    labs(x = "Hygiene (Day 1 of Festival)", y = "Frequency") + 
    theme(legend.position="none")

Pasted image 20241028214955.png

4.7. Boxplots (box-whisker diagrams)

ggplot(festivalData2, aes(x = gender, y = day1, fill = gender)) + 
    geom_boxplot(aes(fill = gender)) +
    scale_fill_manual(values = c("violet", "green"))  +
    labs(x = "Hygiene (Day 1 of Festival)", y = "Frequency") + 
    theme(legend.position="none")

Pasted image 20241028215858.png

4.8. Density plots

4.9. Graphing means

4.9.1. Bar charts and error bars

4.9.1.1. Bar charts for one independent variable

#set cran for r code running in obsidian 
options(repos = c(CRAN = "https://cran.r-project.org"))

#Set the working directory (you will need to edit this to be the directory where you have stored the data files for this Chapter)
setwd("E:/R/dsur/C4exploring_data_with_graphs")

#Initiate packages
library(ggplot2)
library(reshape)
library(plyr)

# import datasets
chickFlick = read.delim("ChickFlick.dat.txt",  header = TRUE)

# if we want means then we have no choice but to dive in stat., here we are using stat_summary to summarize the specified functions
ggplot(chickFlick, aes(x = film, y = arousal)) +
    stat_summary(fun = mean, geom = "bar", fill = "white", color = "black")

Pasted image 20241029115535.png

ggplot(chickFlick, aes(x = film, y = arousal)) +
    stat_summary(fun = mean, geom = "bar", fill = "white", color = "black") +
    stat_summary(fun.data = mean_cl_normal, geom = "pointrange", color = "red")

Pasted image 20241029115551.png

ggplot(chickFlick, aes(x = film, y = arousal)) +
    stat_summary(fun = mean, geom = "bar", fill = "white", color = "black") +
    stat_summary(fun.data = mean_cl_normal, geom = "pointrange", color = "red") +
    xlab("Film") +
    ylab("Mean arousal")

Pasted image 20241029115606.png

4.9.1.2. Bar charts for several independent variables

ggplot(chickFlick, aes(x = film, y = arousal, fill = gender)) +
    stat_summary(fun = mean, geom = "bar", position = "dodge")

Pasted image 20241029120918.png
- 否则，会变成这样
Pasted image 20241029120951.png
- 加上 errorbars, 这里的 position_dodge(width = 0.9) 是 errorbar 的竖线之间的宽度，width = 0.2 是 errorbar 的两端横线的宽度

ggplot(chickFlick, aes(x = film, y = arousal, fill = gender)) +
    stat_summary(fun = mean, geom = "bar", position = "dodge") +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position = position_dodge(width = 0.9), width = 0.2)

Pasted image 20241029121732.png
- 最后加上 labels

ggplot(chickFlick, aes(x = film, y = arousal, fill = gender)) +
    stat_summary(fun = mean, geom = "bar", position = "dodge") +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position = position_dodge(width = 0.9), width = 0.2) +
    xlab("Film") +
    ylab("Mean arousal")

Pasted image 20241029122158.png

ggplot(chickFlick, aes(x = film, y = arousal, fill = film)) +
    stat_summary(fun = mean, geom = "bar", position = "dodge") +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2) +
    facet_wrap(~gender) +
    xlab("Film") +
    ylab("Mean arousal") +
    theme(legend.position = "none")

Pasted image 20241029123004.png

4.9.2. Line graphs

4.9.2.1. Line graphs of a single independent variable

#Initiate ggplot2
library(ggplot2)
library(reshape)
library(plyr)

#Set wd
setwd("E:/R/dsur/C4exploring_data_with_graphs")

#Import dataset
hiccupsData <- read.delim("Hiccups.dat.txt",  header = TRUE)

# transfer the wide format to long format so that ggplot2 can use
hiccups<-stack(hiccupsData)
# rename the column names
names(hiccups)<-c("Hiccups","Intervention")

# transfer "Intervention" as a factor, meanwhile rename it as "Intervention_Factor"
## to avoid duplicated troubles, use unique()
hiccups$Intervention_Factor<-factor(hiccups$Intervention, levels = unique(hiccups$Intervention))

# we can plot the means of each intervention
ggplot(hiccups, aes(x = Intervention_Factor, y = Hiccups)) + 
    stat_summary(fun = mean, geom = "point")
# we want to add a line to connect the means point to show the difference
ggplot(hiccups, aes(x = Intervention_Factor, y = Hiccups)) + 
    stat_summary(fun = mean, geom = "point") +
    stat_summary(fun = mean, geom = "line", aes(group = 1)) #aes(group = 1) means to connect 4 means points as a line, but it doesnt matter about the number, just put a constant in 你可以在 `aes()` 内添加 `group = 1`，这会让 `geom_line()` 将所有数据点视为一个组，从而绘制出正确的连线
# we want to modify the line color and linetype
ggplot(hiccups, aes(x = Intervention_Factor, y = Hiccups)) + 
    stat_summary(fun = mean, geom = "point") +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "blue", linetype = "dashed")
# we want to add errorbars
ggplot(hiccups, aes(x = Intervention_Factor, y = Hiccups)) + 
    stat_summary(fun = mean, geom = "point") +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "blue", linetype = "dashed") +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "green", width = 0.2)
# we want to add xlab and y lab
ggplot(hiccups, aes(x = Intervention_Factor, y = Hiccups)) + 
    stat_summary(fun = mean, geom = "point") +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "blue", linetype = "dashed") +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "green", width = 0.2) +
    xlab("Intervention") +
    ylab("Mean number of hiccups")

Pasted image 20241030085041.png

hiccups$Intervention_Factor<-factor(hiccups$Intervention, levels(hiccups$Intervention)[c(1, 4, 2, 3)])

# we can see the order of representation has been changed
ggplot(hiccups, aes(x = Intervention_Factor, y = Hiccups)) + 
    stat_summary(fun = mean, geom = "point") +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "blue", linetype = "dashed") +
    stat_summary(fun.data = mean_cl_normal, geom = "errorbar", color = "green", width = 0.2) +
    xlab("Intervention") +
    ylab("Mean number of hiccups")

Pasted image 20241030085440.png

4.9.2.2. Line graphs for several independent variables

#Initiate ggplot2
library(ggplot2)
library(reshape)
library(plyr)

#Set wd
setwd("E:/R/dsur/C4exploring_data_with_graphs")

#Import dataset
textData <- read.delim("TextMessages.dat.txt",  header = TRUE)
head(textData)

textData$id <- row(textData[1])

textData$Group_factor <- factor(textData$Group, levels = unique(textData$Group))

textMessages <- reshape(textData, # 要转换的原始数据框
                        idvar = c ("id", "Group_factor"),  # 用于唯一标识每个参与者的列（id和Group）
                        varying = c("Baseline", "Six_months"), # 需要转换的宽格式列，这些列将合并为一列
                        v.names = "Grammar_Score", # 新生成的列名，用于存储合并后的“语法分数”
                        timevar = "Time", # 新生成的时间变量列名，用于区分数据时间点
                        times = c(0:1), # 指定时间点，0表示Baseline，1表示Six_months；这里的数字不重要，若写成 times = c(1:2)那生成的数据框就是用1表示Baseline, 2表示Six_months
                        direction = "long") # 指定转换方向为长格式
### transfer Time as a factor
textMessages$Time_factor <- factor(textMessages$Time, labels = c("Baseline", "After"))

# textMessages1 <- reshape(textData, 
#                         idvar = c ("id", "Group_factor"),
#                         varying = c("Baseline", "Six_months"),
#                         v.names = "Grammar_Score", 
#                         timevar = "Time", 
#                         times = c("Baseline", "After"),
#                         direction = "long")
# textMessages1$Time_factor <- factor(textMessages1$Time, levels = unique(textMessages1$Time))

# textMessages <- melt(textData, id = c("id", "Group"), measured = c("Baseline", "Six_months"))
# names(textMessages) <- c("id", "Group", "Time", "Grammar_Score")
# textMessages$Time <- factor(textMessages$Time, labels = c("Baseline", "After"))

#plot the line graphs
ggplot(textMessages, aes(x = Time_factor, y = Grammar_Score, color = Group_factor)) +
    stat_summary(fun = mean, geom = "point", shape = 9, size = 2) +
    stat_summary(fun = mean, geom = "line", aes(group = Group_factor, linetype = Group_factor)) + # add line to connect points, 不同的组用不同的连接线样式
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.5, size = 1) + # add errorbars using fun.data = mean_cl_boot
    xlab("Time") +
    ylab("Mean Grammar Score")

4.10. Themes and options

ggplot(textMessages, aes(x = Time_factor, y = Grammar_Score, color = Group_factor)) +
    stat_summary(fun = mean, geom = "point", shape = 9, size = 2) +
    stat_summary(fun = mean, geom = "line", aes(group = Group_factor, linetype = Group_factor)) + # add line to connect points, 不同的组用不同的连接线样式
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.5, size = 1) + # add errorbars using fun.data = mean_cl_boot
    xlab("Time") +
    ylab("Mean Grammar Score") +
    theme(panel.grid.major = element_line(color = "violet", linetype = 2)) + #背景的大格子线条颜色和样式
    theme(panel.grid.minor = element_line(color = "purple", linetype = 3)) + #背景的小格子线条颜色和样式
    theme(panel.background = element_rect(fill = "pink")) + #背景填充颜色
    theme(axis.line = element_line(color = "black", linetype = 2, size = 1.5)) #坐标轴线条的颜色和样式

Pasted image 20241031120813.png

Smart Alex's tasks

Task 1

lecturerData = read.delim("Lecturer Data.dat.txt", header = TRUE)

lecturerData$job_factor <- factor(lecturerData$job, labels = c("lecturer", "student"))

friends_comparison <-
    ggplot(lecturerData, aes(x = job_factor, y = friends, color = job_factor)) +
        stat_summary(fun = mean, geom = "bar") +
        stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.5) +
        xlab("Job")+
        ylab("Number of Friends") +
        theme(legend.position = "none")

Pasted image 20241031171227.png

alcohol_comparison <-
    ggplot(lecturerData, aes(x = job_factor, y = alcohol, color = job_factor)) +
        stat_summary(fun = mean, geom = "bar") +
        stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.5) +
        xlab("Job")+
        ylab("Alcohol Consumption") +
        theme(legend.position = "none")

Pasted image 20241031171306.png

income_comparison <-
    ggplot(lecturerData, aes(x = job_factor, y = income, color = job_factor)) +
        stat_summary(fun = mean, geom = "point", size = 2) +
        stat_summary(fun = mean, geom = "line", aes(group = 2), color = "blue", linetype = 2) +
        stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.5) +
        xlab("Job")+
        ylab("Income") +
        theme(legend.position = "none")

Pasted image 20241031171355.png

neurotic_comparison <-
    ggplot(lecturerData, aes(x = job_factor, y = neurotic, color = job_factor)) +
        stat_summary(fun = mean, geom = "point", size = 2) +
        stat_summary(fun = mean, geom = "line", aes(group = 2), color = "blue", linetype = 2) +
        stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.5) +
        xlab("Job")+
        ylab("Neurotic Measure") +
        theme(legend.position = "none")

Pasted image 20241031171428.png

scatter_alcohol_neurotic <-
    ggplot(lecturerData, aes(x = neurotic, y = alcohol, group = job_factor, color = job_factor)) +
        geom_point() +
        geom_smooth(method = "lm", aes(fill = job_factor), alpha = 0.1) +
        scale_fill_manual(values = c("#ffb385", "#000055")) + #填充CI的颜色
        scale_color_manual(values = c("#ffb385", "#000055")) +  #填充拟合线的颜色
        xlab("Neurotic Measure") +
        ylab("Alcohol Consumption")

Pasted image 20241031171514.png

Task 2

infidelityData <- read.csv("Infidelity.csv", header = TRUE)

infidelityData$id <- row(infidelityData[1])
## make Gender as a factor
infidelityData$Gender_factor <- factor(infidelityData$Gender, levels = unique(infidelityData$Gender))
infidelity_long <- reshape(infidelityData,
                           idvar = c("id", "Gender_factor"),
                           varying = c("Partner", "Self"),
                           v.names = "Bullet_number",
                           timevar = "Target",
                           times = c(0:1), #"Partner--0" "Self--1"
                           direction = "long")

infidelity_long$Target_factor <- factor(infidelity_long$Target, labels = c("Partner", "Self"))

ggplot(infidelity_long, aes(x = Target_factor, y = Bullet_number, group = Gender_factor, color = Gender_factor)) +
    stat_summary(fun = mean, geom = "bar", position = "dodge", aes(fill = Gender_factor)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.9), width = 0.2) +
    scale_fill_manual(values = c("#ffb385", "#000055")) + #填充色
    #scale_color_manual(values = c("red", "green")) +  #轮廓线颜色
    xlab("Targets") +
    ylab("Number of Bullets")

Pasted image 20241031172258.png

C5. Exploring assumptions

5.1. What will this chapter tell me?

5.3. Assumptions of parametric data

5.4. Packages used in this chapter

install.packages("car")
install.packages("ggplot2")
install.packages("pastecs")
install.packages("psych")

library(car)
library(ggplot2)
library(pastecs)
library(psych)

5.5. The assumption of normality

5.5.1. Checking normality visually

setwd("E:/R/dsur/C5exploring_assumptions")
imageDirectory <- "E:/R/dsur/C5exploring_assumptions/images"

#Read in the download data:

dlf <- read.delim("DownloadFestival.dat.txt", header=TRUE)
names(dlf)
View(dlf) 
#day1列有一个20 .02的异常值 outlier，我们要把它给改成 NA
#Remove the outlier from the day1 hygiene score
dlf$day1 <- ifelse(dlf$day1 > 20, NA, dlf$day1)

hist.day11 <- ggplot(dlf, aes(day1)) +
    geom_histogram(color = "black", fill = "white") +
    theme(legend.position = "none") +
    stat_function(fun = dnorm, args = list(mean = mean(dlf$day1, na.rm = TRUE), sd = sd(dlf$day1, na.rm = TRUE)), color = "blue", linewidth = 0.5) +
    labs(x = "Hygiene score on day 1", y = "Count")
hist.day11

Pasted image 20250107121441.png

hist.day1 <- ggplot(dlf, aes(day1)) + 
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    theme(legend.position = "none") +
    labs(x = "Hygiene score on day 1", y = "Density")
hist.day1

Pasted image 20250107121205.png

#Histogram for day 2:
hist.day2 <- ggplot(dlf, aes(day2)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    theme(legend.position = "none") +
    labs(x = "Hygiene score on day 2", y = "Density")
hist.day2

#Histogram for day 3:
hist.day3 <- ggplot(dlf, aes(day3)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    theme(legend.position = "none") +
    labs(x = "Hygiene score on day 2", y = "Density")
hist.day3

#Histogram for day1 with a norm curve:
hist.day1 <- ggplot(dlf, aes(day1)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(dlf$day1, na.rm = TRUE), sd = sd(dlf$day1, na.rm = TRUE)), color = "blue", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "Hygiene score on day 1", y = "Density")
hist.day1
ggsave(file = paste(imageDirectory,"05_dlf_day1_hist.png",sep="/"))

Pasted image 20250107121651.png

#Histogram for day2 with a norm curve:
hist.day2 <- ggplot(dlf, aes(day2)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(dlf$day2, na.rm = TRUE), sd = sd(dlf$day2, na.rm = TRUE)), color = "blue", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "Hygiene score on day 2", y = "Density")
hist.day2

Pasted image 20250107121733.png

#Histogram for day3 with a norm curve:
hist.day3 <- ggplot(dlf, aes(day3)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(dlf$day3, na.rm = TRUE), sd = sd(dlf$day3, na.rm = TRUE)), color = "blue", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "Hygiene score on day 3", y = "Density")
hist.day3

Pasted image 20250107121849.png

#Q-Q plot for day 1:
qqnorm(dlf$day1)
qqline(dlf$day1, col = "red")

Pasted image 20250107122005.png

#Q-Q plot for day 2:
qqnorm(dlf$day2)
qqline(dlf$day2, col = "red")

Pasted image 20250107122035.png

#Q-Q plot of the hygiene scores on day 3:
qqnorm(dlf$day3)
qqline(dlf$day3, col = "red")

Pasted image 20250107122103.png

5.5.2. Quantify normality with numbers

zskewness=S−0SEskewnesszkurtosis=K−0SEkurtosis

- 若绝对值>1.96, 则在 p<0.05 的情况下 sign.
- 若绝对值>2.58, 则在 p<0.01 的情况下 sign.
- 若绝对值>3.29, 则在 p<0.001 的情况下 sign.
- 样本量较小的情况，用 p=0.05 就够了；但是当 sample size≥200 时，用 p=0.01 甚至 p=0.001 都不为过
- 【注意：当 sample size ≥ 200 时，先画图出来看 shape of the distribution 很重要，而不是一上来就直接计算 sign.】

### 5.5.2. Quantifying normality with numbers
library(psych)		#load the psych library, if you haven't already, for the describe() function.

#Using the describe() function for a single variable.
describe(dlf$day1)

#Two alternative ways to describe multiple variables.
describe(cbind(dlf$day1, dlf$day2, dlf$day3)) # cbind()把需要的数据列按列拼接，如果直接写describe(dlf)也可以得到想要的结果，但是会把gender和ticknumb列也描述出来
describe(dlf[,c("day1", "day2", "day3")])

#Using the stat.desc() function for a single variable
library(pastecs) #stat.desc() function is from the pastecs package
stat.desc(dlf$day1, basic = FALSE, norm = TRUE)

#Two alternative ways to describe multiple variables.
stat.desc(cbind(dlf$day1, dlf$day2, dlf$day3), basic = FALSE, norm = TRUE)
stat.desc(dlf[, c("day1", "day2", "day3")], basic = FALSE, norm = TRUE)

round(stat.desc(dlf[, c("day1", "day2", "day3")], basic = FALSE, norm = TRUE), digits = 3)

5.5.3. Exploring groups of data

5.5.3.1. Running the analysis for all data

#Read in R exam data.
rexam <- read.delim("rexam.dat.txt", header=TRUE)

#Set the variable uni to be a factor: 0 means Duncetown Uni, 1 means Sussex Uni
rexam$uni<-factor(rexam$uni, levels = c(0:1), labels = c("Duncetown University", "Sussex University"))

#Self-test task
#descriptive stats of the rexam dataset
describe(rexam[, c("exam", "computer", "numeracy", "lectures")])
stat.desc(cbind(rexam$exam, rexam$computer, rexam$numeracy, rexam$lectures), basic = FALSE, norm = TRUE)
round(stat.desc(rexam[, c("exam", "computer", "numeracy", "lectures")], basic = FALSE, norm = TRUE), digits = 3)

round(stat.desc(rexam[, c("exam", "computer", "numeracy", "lectures")], basic = FALSE, norm = TRUE), digits = 3)
exam computer numeracy lectures
median 60.000 51.500 4.000 62.000
mean 58.100 50.710 4.850 59.765
SE.mean 2.132 0.826 0.271 2.168
CI.mean.0.95 4.229 1.639 0.537 4.303
var 454.354 68.228 7.321 470.230
std.dev 21.316 8.260 2.706 21.685
coef.var 0.367 0.163 0.558 0.363
skewness -0.104 -0.169 0.933 -0.410
skew.2SE -0.215 -0.350 1.932 -0.849
kurtosis -1.148 0.221 0.763 -0.285
kurt.2SE -1.200 0.231 0.798 -0.298
normtest.W 0.961 0.987 0.924 0.977
normtest.p 0.005 0.441 0.000 0.077

然后，我们画直方图

#draw histograms
hist.exam <- ggplot(rexam, aes(exam)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(rexam$exam, na.rm = TRUE), sd = sd(rexam$exam, na.rm = TRUE)), color = "#fc00fc", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "First year exam scores", y = "Density")
hist.exam

Pasted image 20250220160851.png

hist.computer <- ggplot(rexam, aes(computer)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(rexam$computer, na.rm = TRUE), sd = sd(rexam$computer, na.rm = TRUE)), color = "#fc00fc", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "Computer Literacy", y = "Density")
hist.computer

Pasted image 20250220160929.png

hist.numeracy <- ggplot(rexam, aes(numeracy)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(rexam$numeracy, na.rm = TRUE), sd = sd(rexam$numeracy, na.rm = TRUE)), color = "#fc00fc", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "Numerical ability", y = "Density")
hist.numeracy

Pasted image 20250220161002.png

hist.lectures <- ggplot(rexam, aes(lectures)) +
    geom_histogram(aes(y = after_stat(density)), color = "black", fill = "white") + 
    stat_function(fun = dnorm, args = list(mean = mean(rexam$lectures, na.rm = TRUE), sd = sd(rexam$lectures, na.rm = TRUE)), color = "#fc00fc", linewidth = 0.5) +
    theme(legend.position = "none") +
    labs(x = "Lectures attended", y = "Density")
hist.lectures

Pasted image 20250220161029.png

5.5.3.2. Running the analysis for different groups

#若想得到不同组别的多个变量的descriptives，可以用cbind()来指定多个变量
by(cbind(data=rexam$exam, data=rexam$numeracy), rexam$uni, describe)
by(rexam[, c("exam", "numeracy")], rexam$uni, stat.desc, basic = FALSE, norm = TRUE)
#可以看到Sussex Uni的学生的exam平均成绩比Dunceton Uni的高36%，numeracy的分数也比较高

#Use describe for four variables in the rexam dataframe.
describe(cbind(rexam$exam, rexam$computer, rexam$lectures, rexam$numeracy))

#Use by() to get descriptives for four variables, split by uni
by(data = cbind(rexam$exam, rexam$computer, rexam$lectures, rexam$numeracy), rexam$uni, describe)
by(rexam[, c("exam", "computer", "lectures", "numeracy")], rexam$uni, stat.desc, basic = FALSE, norm = TRUE)
## round the results with 3 digits
by(
    data = rexam[, c("exam", "computer", "lectures", "numeracy")], 
    INDICES = rexam$uni, 
    FUN = function(subset) {
        round(
            stat.desc(subset, basic = FALSE, norm = TRUE),
            digits = 3
        )
    }
)

#Self test:
#Use by() to get descriptives for computer literacy and percentage of lectures attended, split by uni
by(cbind(data=rexam$computer, data=rexam$lectures), rexam$uni, describe)
by(rexam[, c("computer", "lectures")], rexam$uni, stat.desc, basic = FALSE, norm = TRUE)

#using subset to plot histograms for different groups:
dunceData<-subset(rexam, rexam$uni=="Duncetown University")
sussexData<-subset(rexam, rexam$uni=="Sussex University")

hist.numeracy.duncetown <- 
    ggplot(dunceData, aes(numeracy)) + 
    theme(legend.position = "none") + 
    geom_histogram(aes(y = ..density..), fill = "white", colour = "black", binwidth = 1) + 
    labs(x = "Numeracy Score", y = "Density") + 
    stat_function(fun=dnorm, args=list(mean = mean(dunceData$numeracy, na.rm = TRUE), sd = sd(dunceData$numeracy, na.rm = TRUE)), colour = "red", linewidth = 1)
hist.numeracy.duncetown

Pasted image 20250220162416.png

hist.exam.duncetown <- 
    ggplot(dunceData, aes(exam)) + 
    theme(legend.position = "none") + 
    geom_histogram(aes(y = ..density..), fill = "white", colour = "black") + 
    labs(x = "First Year Exam Score", y = "Density") + 
    stat_function(fun=dnorm, args=list(mean = mean(dunceData$exam, na.rm = TRUE), sd = sd(dunceData$exam, na.rm = TRUE)), colour = "red", linewidth = 1)
hist.exam.duncetown

Pasted image 20250220162447.png

hist.numeracy.sussex <- 
    ggplot(sussexData, aes(numeracy)) + 
    theme(legend.position = "none") + 
    geom_histogram(aes(y = ..density..), fill = "white", colour = "black", binwidth = 1) + 
    labs(x = "Numeracy Score", y = "Density") + 
    stat_function(fun=dnorm, args=list(mean = mean(sussexData $numeracy, na.rm = TRUE), sd = sd(sussexData$numeracy, na.rm = TRUE)), colour = "red", linewidth = 1)
hist.numeracy.sussex

Pasted image 20250220162517.png

hist.exam.sussex <- 
    ggplot(sussexData, aes(exam)) + 
    theme(legend.position = "none") + 
    geom_histogram(aes(y = ..density..), fill = "white", colour = "black") + 
    labs(x = "First Year Exam Score", y = "Density") + 
    stat_function(fun=dnorm, args=list(mean = mean(sussexData$exam, na.rm = TRUE), sd = sd(sussexData$exam, na.rm = TRUE)), colour = "red", linewidth = 1)
hist.exam.sussex

Pasted image 20250220162554.png

5.6. Testing whether a distribution is normal P217

5.6.1. Doing the Shapiro-Wilk test in R

N.B.: Not simply believe the test results, you should also plot the data (in histograms or qqplots) to see if it is normal

#Shapiro-Wilks test for exam and numeracy for whole sample
shapiro.test(rexam$exam)
shapiro.test(rexam$numeracy)
#结果显示，exam和numeracy的shapiro-wilk test结果的p-value都显著< 0.05,说明这两个数据集都不是正态分布的，这也反映了exam 数据有两个modes，numeracy数据是positive skewed的事实

#但是，我们要是按组别来分别对不同uni的数据进行shapiro-wilk test，结果可能就与一饼巴的不同了
#Shapiro-Wilks test for exam and numeracy split by university
by(rexam$exam, rexam$uni, shapiro.test)
#可以看到这时两个uni的exam的test结果都> 0.05 说明数据分布都是normal的；这一点在组间比较中很重要

by(rexam$numeracy, rexam$uni, shapiro.test)
#而两个uni的numeracy数据却non-normal

#qqplots for the two variables
qqPlot(rexam$exam)
ggsave(file = paste(imageDirectory,"05 exam QQ.png",sep="/"))
qqPlot(rexam$numeracy)
ggsave(file = paste(imageDirectory,"05 numeracy QQ.png",sep="/"))
#从这两个qqplot我们能看到总体数据的exam和numeracy的数据分布是non-normal的，因为这些数据点没有与线重合

qqPlot(dunceData$exam)
qqPlot(sussexData$exam)
#而我们单看某个uni的exam数据，就比较normal了，所以要检查数据是否normal，要结合图与test result 一起看

5.6.2. Reporting the Shapiro-Wilk test

shapiro.test(rexam$exam)
#可以这样报告test结果：The percentage on the R exam, W = 0.96, P = .005, was significantly non-normal

Here2day

C6.

C18.