20211122_In class_SCL
Shang Chi Lee
2021/11/22
In class1
how to implement ‘small multiples’ in base graphics given the 4 different diets of the ChickWeight{datasets} example. Adapt the script to produce a plot of 5 panels in which each panel shows a histogram of IQ for each of 5 classes with over 30 pupils in the nlschools{MASS} dataset.
dta <- ChickWeight
dta_diet <- split(dta, dta$Diet)
par(mfrow=c(2, 2), mar=c(4, 4, 0, 0))
lapply(dta_diet, function(x) {
plot(x$weight ~ x$Time,
xlab="Time (day)",
ylab="Weight (gm)")
legend('topleft',
paste("Diet", x$Diet[1], sep="="),
bty='n')}
)
## $`1`
## $`1`$rect
## $`1`$rect$w
## [1] 5.821281
##
## $`1`$rect$h
## [1] 52.72942
##
## $`1`$rect$left
## [1] -0.84
##
## $`1`$rect$top
## [1] 315.8
##
##
## $`1`$text
## $`1`$text$x
## [1] 1.151298
##
## $`1`$text$y
## [1] 289.4353
##
##
##
## $`2`
## $`2`$rect
## $`2`$rect$w
## [1] 5.821281
##
## $`2`$rect$h
## [1] 57.02589
##
## $`2`$rect$left
## [1] -0.84
##
## $`2`$rect$top
## [1] 342.68
##
##
## $`2`$text
## $`2`$text$x
## [1] 1.151298
##
## $`2`$text$y
## [1] 314.1671
##
##
##
## $`3`
## $`3`$rect
## $`3`$rect$w
## [1] 5.821281
##
## $`3`$rect$h
## [1] 65.22824
##
## $`3`$rect$left
## [1] -0.84
##
## $`3`$rect$top
## [1] 386.36
##
##
## $`3`$text
## $`3`$text$x
## [1] 1.151298
##
## $`3`$text$y
## [1] 353.7459
##
##
##
## $`4`
## $`4`$rect
## $`4`$rect$w
## [1] 5.821281
##
## $`4`$rect$h
## [1] 55.26824
##
## $`4`$rect$left
## [1] -0.84
##
## $`4`$rect$top
## [1] 333.32
##
##
## $`4`$text
## $`4`$text$x
## [1] 1.151298
##
## $`4`$text$y
## [1] 305.6859summary(dta)## weight Time Chick Diet
## Min. : 35.0 Min. : 0.00 13 : 12 1:220
## 1st Qu.: 63.0 1st Qu.: 4.00 9 : 12 2:120
## Median :103.0 Median :10.00 20 : 12 3:120
## Mean :121.8 Mean :10.72 10 : 12 4:118
## 3rd Qu.:163.8 3rd Qu.:16.00 17 : 12
## Max. :373.0 Max. :21.00 19 : 12
## (Other):506head(dta )## weight Time Chick Diet
## 1 42 0 1 1
## 2 51 2 1 1
## 3 59 4 1 1
## 4 64 6 1 1
## 5 76 8 1 1
## 6 93 10 1 1data(nlschools, package="MASS")
data1<-nlschools#顯示資料類型和結構
str(data1)## 'data.frame': 2287 obs. of 6 variables:
## $ lang : int 46 45 33 46 20 30 30 57 36 36 ...
## $ IQ : num 15 14.5 9.5 11 8 9.5 9.5 13 9.5 11 ...
## $ class: Factor w/ 133 levels "180","280","1082",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ GS : int 29 29 29 29 29 29 29 29 29 29 ...
## $ SES : int 23 10 15 23 10 10 23 10 13 15 ...
## $ COMB : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...head(data1 )## lang IQ class GS SES COMB
## 1 46 15.0 180 29 23 0
## 2 45 14.5 180 29 10 0
## 3 33 9.5 180 29 15 0
## 4 46 11.0 180 29 23 0
## 5 20 8.0 180 29 10 0
## 6 30 9.5 180 29 10 0#確認描述統計 有6個班級加上1個其他，共7個
summary(data1)## lang IQ class GS SES
## Min. : 9.00 Min. : 4.00 15580 : 33 Min. :10.00 Min. :10.00
## 1st Qu.:35.00 1st Qu.:10.50 5480 : 31 1st Qu.:23.00 1st Qu.:20.00
## Median :42.00 Median :12.00 15980 : 31 Median :27.00 Median :27.00
## Mean :40.93 Mean :11.83 16180 : 31 Mean :26.51 Mean :27.81
## 3rd Qu.:48.00 3rd Qu.:13.00 18380 : 31 3rd Qu.:31.00 3rd Qu.:35.00
## Max. :58.00 Max. :18.00 5580 : 30 Max. :39.00 Max. :50.00
## (Other):2100
## COMB
## 0:1658
## 1: 629
##
##
##
##
## pupils30 <- lapply(data1, function(x) nrow(x[])>30)
which(pupils30 == TRUE)## named integer(0)dat1 <- split(data1, data1$class)
par(mfrow=c(1, 5), mar=c(4, 4, 1, 1))
# histogram
lapply(dat1 [c("5480","15580","15980","16180","18380")],
function(x) {
hist(x$IQ,
xlab = "IQ",
ylab = "Number of student",
main = paste("class", x$class[1], sep=" "))
}
)|>
invisible()
In class2
Deaths per 100,000 from male suicides for 5 age groups and 15 countries are given in the table below. The data set is available as suicides2{HSAUR3}.
data(suicides2, package="HSAUR3")
data2<-suicides2#顯示資料類型和結構
str(data2)## 'data.frame': 15 obs. of 5 variables:
## $ A25.34: num 22 9 22 29 16 28 48 7 8 26 ...
## $ A35.44: num 27 19 19 40 25 35 65 8 11 29 ...
## $ A45.54: num 31 10 21 52 36 41 84 11 18 36 ...
## $ A55.64: num 34 14 31 53 47 49 81 18 20 32 ...
## $ A65.74: num 24 27 49 69 56 52 107 27 28 28 ...head(data2)## A25.34 A35.44 A45.54 A55.64 A65.74
## Canada 22 27 31 34 24
## Israel 9 19 10 14 27
## Japan 22 19 21 31 49
## Austria 29 40 52 53 69
## France 16 25 36 47 56
## Germany 28 35 41 49 52data2<- na.omit(data2[c(1:5)])data2 |> as.data.frame() |> head()## A25.34 A35.44 A45.54 A55.64 A65.74
## Canada 22 27 31 34 24
## Israel 9 19 10 14 27
## Japan 22 19 21 31 49
## Austria 29 40 52 53 69
## France 16 25 36 47 56
## Germany 28 35 41 49 52names(data2)[c(1:5)] <- c("Age25-34", "Age35-44", "Age45-54", "Age55-64" , "Age65-74")boxplot(data2)
#資料轉換(t函數)
dat2<-t(data2)dta2<-as.data.frame.table(data2) head(dta2)## Var1 Var2 Freq.Age25.34 Freq.Age35.44 Freq.Age45.54 Freq.Age55.64
## 1 Canada Age25-34 22 27 31 34
## 2 Israel Age25-34 9 19 10 14
## 3 Japan Age25-34 22 19 21 31
## 4 Austria Age25-34 29 40 52 53
## 5 France Age25-34 16 25 36 47
## 6 Germany Age25-34 28 35 41 49
## Freq.Age65.74
## 1 24
## 2 27
## 3 49
## 4 69
## 5 56
## 6 52names(dta2)<- c("Country","Age_group","Age25-34", "Age35-44", "Age45-54", "Age55-64","Age65-74")head(dta2)## Country Age_group Age25-34 Age35-44 Age45-54 Age55-64 Age65-74
## 1 Canada Age25-34 22 27 31 34 24
## 2 Israel Age25-34 9 19 10 14 27
## 3 Japan Age25-34 22 19 21 31 49
## 4 Austria Age25-34 29 40 52 53 69
## 5 France Age25-34 16 25 36 47 56
## 6 Germany Age25-34 28 35 41 49 52boxplot(dta2[,3:7],
data = dta2,
horizontal=T,
varwidth = T,
frame = FALSE,
cex.axis = 0.6,
xlab="Deaths per 100,000 from male suicides")
abline(v=seq(0, 100, 10), lty=3, col='gray') # grid
stripchart(dta2[,3:7],
data = dta2,
add=T,
col='darkgray', pch=1,
cex=.8,
method='jitter')## Warning in plot.xy(xy.coords(x, y), type = type, ...): "data" 不是一個繪圖參數
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "data" 不是一個繪圖參數
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "data" 不是一個繪圖參數
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "data" 不是一個繪圖參數
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "data" 不是一個繪圖參數
In class3
Use the dataset to replicate the plot
dat3 <- read.table("sat_gpa.txt", header = T)#顯示資料類型和結構
str(dat3)## 'data.frame': 6 obs. of 5 variables:
## $ College: chr "Barnard" "Northwestern" "Bowdoin" "Colby" ...
## $ SAT_No : int 1210 1243 1200 1220 1237 1233
## $ GPA_No : num 3.08 3.1 2.85 2.9 2.7 2.62
## $ SAT_Yes: int 1317 1333 1312 1280 1308 1287
## $ GPA_Yes: num 3.3 3.24 3.12 3.04 2.94 2.8head(dat3)## College SAT_No GPA_No SAT_Yes GPA_Yes
## 1 Barnard 1210 3.08 1317 3.30
## 2 Northwestern 1243 3.10 1333 3.24
## 3 Bowdoin 1200 2.85 1312 3.12
## 4 Colby 1220 2.90 1280 3.04
## 5 Carnegie Mellon 1237 2.70 1308 2.94
## 6 Georgia Tech 1233 2.62 1287 2.80library(devtools)## 載入需要的套件：usethislibrary(ggplot2)#想試著用ggplot做，還未試出，如果用兩張圖重疊再一起?
p <- ggplot(data = dat3,
mapping = aes(
x = SAT_No,
y = College,
label = College))
p + geom_text()
#以下再嘗試 dat3 <- data.frame( select(SAT_No, SAT_Yes, College) %>% filter(SAT %in% c(SAT_No, SAT_Yes)) -> dat31 p <- ggplot(data = dat31, mapping = aes( x = SAT, y = College)) p + geom_line(aes(group = College)) + geom_point(size = 2.0) + geom_text_repel(data = filter(dat31), mapping = aes( x =SAT, y = College, label = College), direction = “x”)
data(organdata, package = “socviz”)
head(organdata)
#參考昶濬作法
par(mar=c(10, 2, 2, 2))
with(dat3, plot(SAT_No ~ GPA_No,
type="n",
xlim=c(1150, 1400),
ylim=c(2.6, 3.4),
xlab="SAT(V+M)",
ylab="First Year GPA"))
# add pair-line
with(dat3[-3,], segments(SAT_No, GPA_No, SAT_Yes, GPA_Yes, lwd=2, col="black"))
segments(x0=1200, y0=2.85, x1=1312, y1=3.12, lwd=4, col=rgb(0,0,0,.6))
# Point
with(dat3, points(x=SAT_Yes, y=GPA_Yes, pch=21, bg="white", cex=2))
with(dat3, points(x=SAT_No, y=GPA_No, pch=21, bg="black", cex=2))
# Text label of College
with(dat3[-3,], text(SAT_Yes+25, GPA_Yes, labels=College, cex=.8))
## font = 1 is the default text font (neither italic nor bold), font = 2 specifies bold face text, font = 3 specifies italic text, and font = 4 specifies both bold and italic text.
with(dat3, text(1312+25, 3.12, labels="Bowdoin", cex=.8, font=2))
# add legend
legend("topleft", c("Submitted SAT Scores", "Did NOT Submit SAT Scores"),
pch=c(21, 19),
bg=c("white","black"),
pt.cex=1.5, bty="n")
# set axis and scale but without label
axis(1, seq(1150, 1400, length=16), labels = F)
axis(2, seq(2.6, 3.4, length=21), labels = F)
# add annotation
mtext("Figure 1.4 The mean SAT coupled with the mean first-year GPA for the class \n of 1999 at six schools shown for those who submitted SAT scores for admis- \n sion and those who did not", side=1, cex=1, adj = 0, line=7)
Exercises1
The data set LSAT{ltm} contains answers by 1,000 persons to 5 test items. Generate a plot similar to the one below:
data(LSAT, package="ltm")
e1<-LSATstr(e1)## 'data.frame': 1000 obs. of 5 variables:
## $ Item 1: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Item 2: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Item 3: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Item 4: int 0 0 0 0 0 0 0 0 0 1 ...
## $ Item 5: int 0 0 0 1 1 1 1 1 1 0 ...head(e1)## Item 1 Item 2 Item 3 Item 4 Item 5
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 1
## 5 0 0 0 0 1
## 6 0 0 0 0 1library(dplyr)##
## 載入套件：'dplyr'## 下列物件被遮斷自 'package:stats':
##
## filter, lag## 下列物件被遮斷自 'package:base':
##
## intersect, setdiff, setequal, unione11 <- e1 |>
dplyr::mutate(Total = rowSums(e1,c(1:5)))dplyr 套件中的基礎函數 mutate() 函數：新增變數 pivot_longer() 長數據 pivot_wider() 橫
e11 <- e1 |>
dplyr::mutate(Total = rowSums(e1,c(1:5)))|>
tidyr::pivot_longer(cols = starts_with("Item"),
names_to = "Item",
values_to = "Score")#計算頻率
e11frq <- as.data.frame(table(e11))head(e11frq)## Total Item Score Freq
## 1 0 Item 1 0 3
## 2 1 Item 1 0 10
## 3 2 Item 1 0 23
## 4 3 Item 1 0 25
## 5 4 Item 1 0 15
## 6 5 Item 1 0 0e11frq1 <- e11frq |>
dplyr::group_by(Total, Item)|>#資料依Total, Item排列
dplyr::mutate(Procorrect=Freq/sum(Freq))|> #計算百分比
subset(Total!=c(0,5)& Score==1)
e11frq1$Total<-as.numeric(as.character(e11frq1$Total))e11frq1## # A tibble: 20 x 5
## # Groups: Total, Item [20]
## Total Item Score Freq Procorrect
## <dbl> <fct> <fct> <int> <dbl>
## 1 1 Item 1 1 10 0.5
## 2 2 Item 1 1 62 0.729
## 3 3 Item 1 1 212 0.895
## 4 4 Item 1 1 342 0.958
## 5 1 Item 2 1 1 0.05
## 6 2 Item 2 1 24 0.282
## 7 3 Item 2 1 109 0.460
## 8 4 Item 2 1 277 0.776
## 9 1 Item 3 1 1 0.05
## 10 2 Item 3 1 7 0.0824
## 11 3 Item 3 1 63 0.266
## 12 4 Item 3 1 184 0.515
## 13 1 Item 4 1 2 0.1
## 14 2 Item 4 1 28 0.329
## 15 3 Item 4 1 139 0.586
## 16 4 Item 4 1 296 0.829
## 17 1 Item 5 1 6 0.3
## 18 2 Item 5 1 49 0.576
## 19 3 Item 5 1 188 0.793
## 20 4 Item 5 1 329 0.922plot(0, 0, bty = "n", type = "n",
xlim = c(0, 5),
ylim = c(0, 1),
xlab = "Total score",
ylab = "Proportional Correct",
)
# point and line
lines(e11frq1$Total[e11frq1$Item == "Item 1"],
e11frq1$Procorrect[e11frq1$Item == "Item 1"],
type = "b",
col=1,
pch = 1)
lines(e11frq1$Total[e11frq1$Item == "Item 2"],
e11frq1$Procorrect[e11frq1$Item == "Item 2"],
type = "b",
col = 2,
pch = 2)
lines(e11frq1$Total[e11frq1$Item == "Item 3"],
e11frq1$Procorrect[e11frq1$Item == "Item 3"],
type = "b",
col = 3,
pch = 3)
lines(e11frq1$Total[e11frq1$Item == "Item 4"],
e11frq1$Procorrect[e11frq1$Item == "Item 4"],
type = "b",
col = 4,
pch = 4)
lines(e11frq1$Total[e11frq1$Item == "Item 5"],
e11frq1$Procorrect[e11frq1$Item == "Item 5"],
type = "b",
col = 5,
pch = 5)
grid()
legend("topleft",
c("Item 1", "Item 2", "Item 3", "Item 4", "Item 5"),
pch = c(1, 2, 3, 4, 5),
col = c(1, 2, 3, 4, 5),
bty = "n", lty = 1)
Exercises2
Complete the R script to create the following plot
# generate weights
wt <- seq(40, 160, length.out=600)
# generate heights
ht <- seq(1.40, 2.00, length.out=600)
# cross wt by ht
wtht <- expand.grid(wt=wt, ht=ht)
# function to compute bmi from wt and ht
bmi_wh <- function(w, h) {w/(h*h)}head(wtht)## wt ht
## 1 40.00000 1.4
## 2 40.20033 1.4
## 3 40.40067 1.4
## 4 40.60100 1.4
## 5 40.80134 1.4
## 6 41.00167 1.4bmi_wh## function(w, h) {w/(h*h)}str(bmi_wh)## function (w, h)
## - attr(*, "srcref")= 'srcref' int [1:8] 12 11 12 34 11 34 12 12
## ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002443e110># generate data matrix 計算BMI矩陣
bmiwtht <- matrix(bmi_wh(wtht$wt, wtht$ht), length(ht), length(wt))#rgb()函數描述顏色
colorbar<-c(rgb(1, 1, 1, 0), rgb( 1, 1, 1, 0),rgb(.94, .50, .50, 0.3))# draw the contour
contour(wt, ht, bmiwtht, bty='n',
levels = c(18.5, 24.9, 30),
drawlabels=FALSE,
ylab="Height (m)",
xlab="Weight (kg)",
main="BMI categories by height and weight",
)
# add grid lines
grid()
# annotate the bmi categories
text(105, 1.8, "Obese", cex=1, srt=45,)
text(92, 1.8, "Overweight", cex=1, srt=45)
text(75, 1.8, "Normal", cex=1, srt=45)
text(55, 1.8, "Underweight", cex=1, srt=55)
#.filled.contour(x, y, z, levels, col)填充為純色
.filled.contour(wt, ht, bmiwtht,levels = c(18.5, 24.9, 30, 100), col= colorbar)
Exercises3
Use the free recall data to replicate the figure
list.files("C:/Nicole/Rstudio IDE code/R Base graphics", pattern = "fr")## [1] "fr10-2.txt" "fr15-2.txt" "fr20-1.txt" "fr20-2.txt" "fr30-1.txt"
## [6] "fr40-1.txt"fr10_2 <- read.table("fr10-2.txt", sep = "", col.names = c(1:10), fill = T, nrows = 1200)
fr10_2["Group"] = "10-2"
fr15_2 <- read.table("fr15-2.txt", sep = "", col.names = c(1:15), fill = T, nrows = 1200)
fr15_2["Group"] = "15-2"
fr20_1 <- read.table("fr20-1.txt", sep = "", col.names = c(1:20), fill = T, nrows = 1200)
fr20_1["Group"] = "20-1"
fr20_2 <- read.table("fr20-2.txt", sep = "", col.names = c(1:20), fill = T, nrows = 1200)
fr20_2["Group"] = "20-2"
fr30_1 <- read.table("fr30-1.txt", sep = "", col.names = c(1:30), fill = T, nrows = 1200)
fr30_1["Group"] = "30-1"
fr40_1 <- read.table("fr40-1.txt", sep = "", col.names = c(1:40), fill = T, nrows = 1200)
fr40_1["Group"] = "40-1"fr_all<-list(fr10_2, fr15_2, fr20_1, fr20_2, fr30_1, fr40_1)#從寬轉為長，共1200筆
fr_all_1 <- lapply(fr_all,
function(x) {y <- x |>
tidyr::pivot_longer(cols = starts_with("X"),
names_to = "Serial",
values_to = "Item")
y <- as.data.frame(table(y$Item)/1200) |> #計算item平均頻率
dplyr::mutate(Group = x$Group[1])
} ) |>
dplyr::bind_rows()
fr_all_1$Group <- factor(fr_all_1$Group)#整理資料
fr_all_1$Group <- factor(fr_all_1$Group)
# remove "88"和 其他非此資料夾應有的資料，如15-2裡不該有大於15的資料
fr_all_2 <- fr_all_1[!fr_all_1$Var1 == "88", ]
fr_all_2 <- fr_all_2[!(fr_all_2$Var1 == "16" & fr_all_2$Group == "15-2"), ]
fr_all_2 <- fr_all_2[!(fr_all_2$Var1 == "31" & fr_all_2$Group == "30-1"), ]
fr_all_2 <- fr_all_2[!(fr_all_2$Var1 == "41" & fr_all_2$Group == "40-1"), ]
fr_all_2<- fr_all_2[!(fr_all_2$Var1 == "50" & fr_all_2$Group == "40-1"), ]
head(fr_all_2)## Var1 Freq Group
## 1 1 0.7025000 10-2
## 2 2 0.5691667 10-2
## 3 3 0.4750000 10-2
## 4 4 0.4633333 10-2
## 5 5 0.4591667 10-2
## 6 6 0.5675000 10-2#學昶濬
# 設置畫圖區域
## xlim 設大一點，最後一個點就不會超出了
plot(0,0, bty="n", xlim = c(0,42), ylim=c(0,1),
xlab = "SERIA POSITION",
ylab = "PROBABILITY OF RECALL",
xaxs = "i", yaxs ="i") # i, 軸線相連；r,軸線有縫隙
# point and line
points(fr_all_2$Var1[fr_all_2$Group == "10-2"],
fr_all_2$Freq[fr_all_2$Group == "10-2"],
pch = 19)
lines(fr_all_2$Var1[fr_all_2$Group == "10-2"],
fr_all_2$Freq[fr_all_2$Group == "10-2"])
points(fr_all_2$Var1[fr_all_2$Group == "15-2"],
fr_all_2$Freq[fr_all_2$Group == "15-2"],
pch = 1)
lines(fr_all_2$Var1[fr_all_2$Group == "15-2"],
fr_all_2$Freq[fr_all_2$Group == "15-2"])
points(fr_all_2$Var1[fr_all_2$Group == "20-2"],
fr_all_2$Freq[fr_all_2$Group == "20-2"],
pch = 19)
lines(fr_all_2$Var1[fr_all_2$Group == "20-2"],
fr_all_2$Freq[fr_all_2$Group == "20-2"])
points(fr_all_2$Var1[fr_all_2$Group == "20-1"],
fr_all_2$Freq[fr_all_2$Group == "20-1"],
pch = 1)
lines(fr_all_2$Var1[fr_all_2$Group == "20-1"],
fr_all_2$Freq[fr_all_2$Group == "20-1"])
points(fr_all_2$Var1[fr_all_2$Group == "30-1"],
fr_all_2$Freq[fr_all_2$Group == "30-1"],
pch = 1)
lines(fr_all_2$Var1[fr_all_2$Group == "30-1"],
fr_all_2$Freq[fr_all_2$Group == "30-1"])
points(fr_all_2$Var1[fr_all_2$Group == "40-1"],
fr_all_2$Freq[fr_all_2$Group == "40-1"],
pch = 19)
lines(fr_all_2$Var1[fr_all_2$Group == "40-1"],
fr_all_2$Freq[fr_all_2$Group == "40-1"])
# text label
text(2, 0.9, "10-2") #text(x,y,lable)
lines(c(2, 6), c(0.85, 0.6)) #lines(c(x1,x2),c(y1,y2))
text(14, 0.65, "15-2")
lines(c(12, 13), c(0.61, 0.63))
text(13, 0.48, "20-2")
lines(c(12, 13), c(0.45, 0.3))
text(18.5, 0.44, "20-1")
lines(c(16, 17), c(0.35, 0.4))
text(23, 0.6, "30-1")
lines(c(21.5, 21, 26), c(0.55, 0.52, 0.39))
text(35, 0.8, "40-1")
lines(c(33.5, 33, 37.5), c(0.77, 0.75, 0.59))
# axis (圖上的軸在這邊設定)
axis(1,seq(0, 40, 10), labels = FALSE)
axis(2,seq(0, 1, 0.1), labels = FALSE)
Exercises4
(Optional) A dataset, raz2005{rmcorr}, contains two repeated measures, on two occasions (Time), of age (in years) and adjusted volume of cerebellar hemispheres from 72 participants. Use it to replicate the following figure in which participants older than 65 years of age at Time One are plotted on a shifted panel and those whose cerebellar volume increased from Time One to Time Two are shown in blue color.
data(raz2005, package="rmcorr")
raz<-raz2005
head(raz)## Participant Time Age Volume
## 1 1 1 24.04637 168.0261
## 2 1 2 29.28124 157.5085
## 3 2 1 21.02195 164.0444
## 4 2 2 25.91324 152.3567
## 5 3 1 23.92362 150.4820
## 6 3 2 29.57359 146.3979#依時間切割
raz1 <- split(raz, raz$Time)#四個邊的順序是下、左、上、右。
par(mar=c(10, 10, 0.5, 0.5))
with(raz1, plot(0,0,type="n",
xlim=c(20, 90),
ylim=c(100, 170),
xlab="Age(year)",
ylab="Cerebellar volume"))
#用loop畫圖
pchlist<-c(1,20)
for (i in 1:2){
points(raz1[[i]]$Age,
raz1[[i]]$Volume,
pch = pchlist[i])}
arrows(raz1[[1]]$Age,raz1[[1]]$Volume,
raz1[[2]]$Age,raz1[[2]]$Volume,
length=0.1, col = "grey")