208329: Lab material
กมลรัตน์ สุภาวรรณ
Kamonrat Suphawan
Department of Statistics
Faculty of Science, Chiang Mai University
Kamonrat Suphawan
Department of Statistics
Faculty of Science, Chiang Mai University
Data files
ปฎิบัติการ: Introduction to R
1. Installing the R and RStudio
R is a programming language for statistical computing and graphics supported by the R Core Team and the R Foundation for Statistical Computing. R is used among data miners, bioinformaticians and statisticians for data analysis and developing statistical software. Users have created packages to augment the functions of the R language.
To work with R, you need to download the program and install it. I recommend that you also install R-Studio. R-Studio is a separate program that makes R easier to use .
Windows
Install R
Download the latest R installer (.exe) for Windows. Install the downloaded file as any other windows app.
Install RStudio
Now that R is installed, you need to download and install RStudio. First download the installer for Windows. Run the installer (.exe file) and follow the instructions.
Mac
Install R
First download the latest release (“R-version.pkg”) of R Save the .pkg file, double-click it to open, and follow the installation instructions. Now that R is installed, you need to download and install RStudio.
Install RStudio
First download the the version for Mac. After downloading, double-click the file to open it, and then drag and drop it to your applications folder.
2. Basic data types
R works with numerous data types. Some of the most basic types to get started are:
- Decimal values like
4.5are called numerics. - Whole numbers like
4are called integers. Integers are also numerics. - Boolean values (
TRUEorFALSE) are called logical. - Text (or string) (e.g. “world”) values are called characters.
Numeric and integer values
# create create variable "a" that is a vector of one number.
a <- 7
a = 7
a## [1] 7# create vectors that consists of multiple numbers.
b <- c(1.25, 2.9, 3.0)
b## [1] 1.25 2.90 3.00# create a regular sequence of whole numbers
c <- 1:10
c## [1] 1 2 3 4 5 6 7 8 9 10# seq(from, to, by) function ro creat sequence of numbers
d <- seq(1,10,0.5)
d## [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0
## [16] 8.5 9.0 9.5 10.0# rep() function to repeat a single number, or a sequence of numbers.
rep(99, times=5)## [1] 99 99 99 99 99rep(c(1,2,3), times=5)## [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3rep(c(1,2,3),each = 5)## [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3Character values
x <- "Monday"
class(x)## [1] "character"y <- c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")
class(y)## [1] "character"Factors
A factor is a nominal (categorical) variable with a set of known
possible values called levels. They can be created using the
as.factor function.
yy <- as.factor(y)
yy## [1] Mon Tue Wed Thu Fri Sat Sun
## Levels: Fri Mon Sat Sun Thu Tue Wed3. Data structure
Matrix
mat0 <- matrix(NA, ncol=3, nrow=2)
mat0## [,1] [,2] [,3]
## [1,] NA NA NA
## [2,] NA NA NAmat1 <- matrix(1:6, ncol=3, nrow=2)
mat1## [,1] [,2] [,3]
## [1,] 1 3 5
## [2,] 2 4 6mat2 <- matrix(1:6, ncol=3, nrow=2,byrow = TRUE)
mat2## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6List
A list is a very flexible container to store data. Each element of a list can contain any type of R object, e.g. a vector, matrix, data.frame, another list, or more complex data types.
lst <- list(d,yy,mat2)
lst## [[1]]
## [1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0
## [16] 8.5 9.0 9.5 10.0
##
## [[2]]
## [1] Mon Tue Wed Thu Fri Sat Sun
## Levels: Fri Mon Sat Sun Thu Tue Wed
##
## [[3]]
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6Data frame
The data.frame is commonly used for statistical data
analysis in R. It is a special type of list that requires that all
elements (variables) have the same length.
# Loading
data(ToothGrowth)
# Print the first 6 rows
head(ToothGrowth,row=6)## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5# Examine data structure of the dataframe
str(ToothGrowth)## 'data.frame': 60 obs. of 3 variables:
## $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
## $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
## $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...4. Read files
To read a file into R, you need to know the full path (directory) name and the name of the file itself.
- Save .xlsx data as .csv file
- Check the directory of data file (right click–> Properties–>Location)
- Use read.csv from base R (Slowest method, but works fine for smaller datasets)
# Set working directory (folder that contains data file)
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
# Read data file and keep it in dt
dt <- read.csv("bank_data.csv", header = TRUE)
head(dt) # show the first six rows of the data## age job marital education balance housing loan duration campaign pdays
## 1 59 admin. married secondary 2343 yes no 1042 1 -1
## 2 56 admin. married secondary 45 no no 1467 1 -1
## 3 41 technician married secondary 1270 yes no 1389 1 -1
## 4 55 services married secondary 2476 yes no 579 1 -1
## 5 54 admin. married tertiary 184 no no 673 2 -1
## 6 42 management single tertiary 0 yes yes 562 2 -1str(dt) # show data structure ## 'data.frame': 1000 obs. of 10 variables:
## $ age : int 59 56 41 55 54 42 56 60 37 28 ...
## $ job : chr "admin." "admin." "technician" "services" ...
## $ marital : chr "married" "married" "married" "married" ...
## $ education: chr "secondary" "secondary" "secondary" "secondary" ...
## $ balance : int 2343 45 1270 2476 184 0 830 545 1 5090 ...
## $ housing : chr "yes" "no" "yes" "yes" ...
## $ loan : chr "no" "no" "no" "no" ...
## $ duration : int 1042 1467 1389 579 673 562 1201 1030 608 1297 ...
## $ campaign : int 1 1 1 1 2 2 1 1 1 3 ...
## $ pdays : int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...Set “chr” variables as “factor”
dt$job <- as.factor(dt$job)
dt$marital <- as.factor(dt$marital)
dt$education <- as.factor(dt$education)
dt$housing <- as.factor(dt$housing)
dt$loan <- as.factor(dt$loan)
#check data structure again
str(dt)5. Numerical presentation for qualitative data or categorical data
summary()function- frequency table with
table()function - contingency table (Two-Way Frequency Tables)
- percentage or proportion with
prop.table()function
# job, marital, education, housing, loan
summary(dt)
#frequency table for one categorical variable
table(dt$job)
table(dt$housing)
prop.table(table(dt$job))
prop.table(table(dt$job))*100
#Two-Way Frequency Tables
table(dt$job,dt$education)
tb <- table(dt$loan,dt$marital)
# Add margin totals to your table.
addmargins(tb)
#Calculate percentage
prop.table(tb)*1006. Numerical presentation for quantitative data
summary()functionmax()andminfunctionmean()functionmedian()function- There is no function to compute mode in R. To compute mode, use the
getmode()function defined below.
#define function to calculate mode
getmode <- function(x) {
u <- unique(x)
tab <- tabulate(match(x, u))
u[tab == max(tab)]
}range()functionquantile(, probs = )function whereprobscan be 0.25, 0.5, 0.75 for first, second, and third quantileIQR()functionsd()andvar()function
#age, balance, duration, campaign, pdays
summary(dt$age)
min(dt$age)
range(dt$age)
mean(dt$age)
getmode(dt$age)
IQR(dt$age)
sd(dt$age)
quantile(dt$age,probs = 0.25)
quantile(dt$age,probs = 0.50)
quantile(dt$age,probs = 0.75)
quantile(dt$age,probs = c(0.25,0.50,0.75))
cv_age <- sd(dt$age)/mean(dt$age)
cv_age
#cv of age by loan (yes,no)
cv1 <- sd(dt$age[dt$loan=="yes"])/mean(dt$age[dt$loan=="yes"])
cv1
cv2 <- sd(dt$age[dt$loan=="no"])/mean(dt$age[dt$loan=="no"])
cv2aggregate()function is used to get the summary statistics of the data by group.
# aggregate: compute summary statistics of data subsets
aggregate(balance~loan, data = dt, FUN=summary)
aggregate(balance~loan+marital, data = dt, FUN=summary)
aggregate(balance~loan+marital, data = dt, FUN=mean)7. Graphical presentation for qualitative data or categorical data
In the case that you have a factor variable, you can display a frequency table or a bar graph.
You can use the good ol’ summary() on a factor variable
to display a frequency table.
Bar graph
freq <- table(dt$job)
freq##
## admin. blue-collar entrepreneur housemaid management
## 96 226 45 25 221
## retired self-employed services student technician
## 44 36 82 8 183
## unemployed unknown
## 31 3barplot(freq)
To obtain a list of a factor variable’s levels, use
levels().
maritallevel <- levels(dt$marital)
maritallevel## [1] "divorced" "married" "single"Stacked Bar Plot with Colors and Legend
counts <- table(dt$marital, dt$loan)
barplot(counts,
legend = rownames(counts),xlab="loan")Grouped Bar Plot
counts <- table(dt$marital, dt$loan)
barplot(counts,col=c("blue","red","green"),
legend = rownames(counts), beside=TRUE,xlab="loan")Pie chart
counts <- table(dt$marital)
pie(counts, main = "Marital pie chart")
8. Graphical presentation for quantitative data
Histogram
- To display a histogram of a vector, use
hist()
hist(dt$age)
Sometimes, you want to change the title of the plot, and also the
label of the horizontal axis and vertical axis. To do that, set the
parameters main, xlab and
ylab
hist(dt$age, main = "Histogram of age", xlab = "age", ylab = "counts")To change the number of bins, set the breaks parameter.
However, R will decide what is the best number of bins close to the
number that your specify. For example, the following code displays a
histogram with the number of bins close to 9:
hist(dt$age, breaks = 20,
main = "Histogram of age",
xlab = "age",
ylab = "counts")Stem and leaf plot
- To display a stem and leaf plot of a vector, use
stem().
stem(dt$duration)##
## The decimal point is 2 digit(s) to the right of the |
##
## 0 | 18890113333334455566666667788899
## 2 | 1333333344555577888901111244456667777788888889999
## 4 | 00000011112222222333333334444445555566666666777777777888888888888899+122
## 6 | 00000000001111111122222222222223333333333333334444444444444444445555+160
## 8 | 00000000001111111111222223333333444444444455555555555556666666666677+103
## 10 | 00000000011111112222222233333333333344444444455555566666666667788888+54
## 12 | 00000111111112233344444444555667778899900011122333344444445556666677
## 14 | 0001122344455556777899900133445556667788
## 16 | 012245667899022234899
## 18 | 1346778885788
## 20 | 2239
## 22 | 7
## 24 | 62
## 26 | 257
## 28 |
## 30 | 98
## 32 | 5
## 34 |
## 36 |
## 38 | 8#compare with its histogram Boxplot
- To display a boxplot of a vector, use
boxplot().
boxplot(dt$age,main="Boxplot of age")
- To display the boxplot horizontally, set the parameter
horizontal = TRUE
boxplot(dt$age, horizontal = TRUE)You can display boxplots of by group variable.
boxplot(dt$age~dt$loan,xlab="loan",ylab="age")boxplot(dt$age~dt$job,xlab="job",ylab="age")Scatterplot
Note that you need two vectors to plot a scatter plot: one for the x-axis, and one for the y-axis.
Suppose that you have two variables: var1 and
var2. To display a scatter plot between var1
and var2, use plot(x = var1, y = var2).
plot(x = dt$age, y = dt$balance,xlab="age",ylab="balance")9. Assignment
Perform the following numerical data presentation for the
bank_data.csv dataset:
- Create frequency table of education, with proportion
- Create two-ways frequency table of job by loan category, with percentage
- Compute summary statistics of balance
- Compute Q1,Q2, and Q3 of balance, also the IQR
- Compute CV of balance by type of housing
- Compute summary statistics of balance by housing
- Compute sd of balance by housing and loan
Perform the following graphical data presentation for the
bank_data.csv dataset:
- How would you present the education data?
- If you would like to examine the number of people with different education and martial status, which graph shall you use to present such information?
- How would you present the balance data?
- If you would like to compare the balance between people with different education, which graph shall you use to present such information?
- If you would like to examine the relationship between age and duration, which graph shall you use to present such information?
ปฎิบัติการ: การวิเคราะห์การถดถอยเชิงเส้นอย่างง่าย (Simple linear regression)
วัตถุประสงค์ นักศึกษาสามารถ
- แสดงการวิเคราะห์การถดถอยเชิงเส้นอย่างได้
- แสดงการจำลองข้อมูลเพื่อตรวจสอบการแจกแจงของตัวประมาณพารามิเตอร์การถดถอยเชิงเส้นได้
Task 1: The Rocket data
Set working directory (Folder ที่มีไฟล์ข้อมูล)ด้วยคำสั่ง setwd()
และ import data ด้วยคำสั่ง read.csv()
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
dat1 <- read.csv("TheRocketData.csv")
#Check data structure
str(dat1)## 'data.frame': 20 obs. of 3 variables:
## $ Observation : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Sher.Strength : num 2159 1678 2316 2061 2208 ...
## $ Age.of.Propellant: num 15.5 23.8 8 17 5.5 ...เพื่อให้สะดวกในการพิมพ์ชื่อตัวแปร กำหนด
x <- dat1$Age.of.Propellant
y <- dat1$Sher.Strengthสร้างแผนภาพการกระจายระหว่าง x และ y
plot(x,y,xlab="Age.of.Propellant (X)",ylab="Sher.Strength(Y)")
คำนวณค่าสัมประสิทธิ์สหสัมพันธ์ของตัวอย่าง และด้วยสัมประสิทธิ์สหสัมพันธ์
จงทำการทดสอบที่ระดับนัยสำคัญ 0.05 ว่า x และ y มีความสัมพันธ์กันหรือไม่ ด้วยคำสั่ง
cor.test()
cor.test(x,y)##
## Pearson's product-moment correlation
##
## data: x and y
## t = -12.86, df = 18, p-value = 1.643e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9802377 -0.8747304
## sample estimates:
## cor
## -0.9496533ทำการวิเคราะห์ข้อมูลดังกล่าวด้ยการวิเคราะห์การถดถอยเชิงเส้นอย่างง่าย ด้วยคำสั่ง
lm()
fit1 <- lm(y ~ x)
summary(fit1)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -215.98 -50.68 28.74 66.61 106.76
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2627.822 44.184 59.48 < 2e-16 ***
## x -37.154 2.889 -12.86 1.64e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 96.11 on 18 degrees of freedom
## Multiple R-squared: 0.9018, Adjusted R-squared: 0.8964
## F-statistic: 165.4 on 1 and 18 DF, p-value: 1.643e-10anova(fit1)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 1527483 1527483 165.38 1.643e-10 ***
## Residuals 18 166255 9236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1หาช่วงความเชื่อมั่นของพารามิเตอร์การถดถอย ด้วยคำสั่ง
confint.lm()
#90% confidence interval
confint.lm(fit1,level=0.90)
#95% confidence interval
confint.lm(fit1,level=0.95)คำนวณค่าส่วนเหลือ (residuals) ด้วยคำสั่ง resid()
res <- resid(fit1)
#residuals diagnostic
par(mfrow=c(2,2))#create 2x2 layout plotting
plot(fit1)
หาค่าพยากรณ์ของตัวแปรตามจากค่าของ x ด้วยคำสั่ง predict()
yhat <- predict(fit1)คำนวณช่วงความเชื่อมั่น 95 % ของค่าเฉลี่ยของ y และค่าพยากรณ์ y สำหรับทุกๆค่าที่เป็นไปได้ของ x
#Confidence intervals for E(Y)
new <- data.frame(x = seq(min(x),max(x),len=100))
pred.cf <- predict(fit1, new, interval = "confidence")
#Prediction intervals for Y
pred.pd <- predict(fit1, new, interval = "prediction") สร้างแผนภาพแสดงข้อมูลรวมถึงค่าพยากรณ์ของ y (\(\hat{y}\)) และช่วงความเชื่อมั่น 95 % ของค่าเฉลี่ยของ y และช่วงพยากรณ์ค่าของ y
plot(x,y,xlab="Age of Propellant",ylab="Sher Strength")
points(x,yhat,col="red") #predicted values
lines(x,yhat,col="red") #fitted regression line
# 95 % confidence interval of E(y) for all possible values of x
lines(new$x,pred.cf[,"lwr"],col="blue",lty=2)
lines(new$x,pred.cf[,"upr"],col="blue",lty=2)
# 95 \% prediction interval of y for all possible values of x
lines(new$x,pred.pd[,"lwr"],col="green",lty=2)
lines(new$x,pred.pd[,"upr"],col="green",lty=2)
Task 2: Simulation study on the distribution of the regression coefficients
Suppose our population of interest consists of 5 pairs of \((X, Y)\) \[\begin{array}{c|ccccc} no. & 1 & 2 & 3 & 4 &5 \\ \hline X & 8 & 12 & 29.5 & 43 & 53 \\ \hline Y & 8.89 & 8.14 & 6.67 & 6.08 &5.85(เปลี่ยนเป็นเลขสามตัวท้ายของรหัสนักศึกษา) \end{array}\]Set the data frame for the population data
popl <- data.frame(X=c(8,12,29.5,43,53),
Y=c(8.89,8.14,6.67,6.08,5.85))
str(popl)## 'data.frame': 5 obs. of 2 variables:
## $ X: num 8 12 29.5 43 53
## $ Y: num 8.89 8.14 6.67 6.08 5.85Fit regression line to the population data. You will obtain the true parameter values \((\beta_{0}, \beta_{1})\). Also, plot the true regression line on the data
fit.pop <- lm(Y ~ X, data=popl)
summary(fit.pop)##
## Call:
## lm(formula = Y ~ X, data = popl)
##
## Residuals:
## 1 2 3 4 5
## 0.3625 -0.1218 -0.4294 -0.1227 0.3115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.058858 0.335722 26.983 0.000112 ***
## X -0.066421 0.009912 -6.701 0.006780 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3842 on 3 degrees of freedom
## Multiple R-squared: 0.9374, Adjusted R-squared: 0.9165
## F-statistic: 44.91 on 1 and 3 DF, p-value: 0.00678plot(popl$X,popl$Y,xlab="X",ylab="Y")
abline(fit.pop)
Randomly sample with size \(n=3\) pairs from the population and store all possible sets. You will obtain 10 possible datasets. Then, fit the linear regression to each dataset.
# Example: set 1
x1 <- c(8,12,29.5)
y1 <- c(8.89,8.14,6.67)
fit1 <- lm(y1~x1)
summary(fit1)##
## Call:
## lm(formula = y1 ~ x1)
##
## Residuals:
## 1 2 3
## 0.16162 -0.19856 0.03694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.50804 0.30325 31.353 0.0203 *
## x1 -0.09746 0.01600 -6.093 0.1036
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2587 on 1 degrees of freedom
## Multiple R-squared: 0.9738, Adjusted R-squared: 0.9475
## F-statistic: 37.12 on 1 and 1 DF, p-value: 0.1036Examine sampling distribution of \(b_{0}\) and \(b_{1}\) and compare to the parameter values.
#TO DOปฎิบัติการ: การวิเคราะห์การถดถอยเมื่อมีตัวแปรหุ่นในตัวแบบ
วัตถุประสงค์ นักศึกษาสามารถ
- แสดงการแปลงตัวแปรอิสระที่เป็นตัวแปรเชิงคุณภาพให้เป็นตัวแปรหุ่นได้
- แสดงการวิเคราะห์การถดถอยในการพยากรณ์ตัวแปรตาม y เมื่อตัวแปรอิสระมีทั้งตัวแปรเชิงคุณภาพและเชิงปริมาณ พร้อมทั้งอธิบายความหมายของสมการที่ได้ได้
- แสดงการวิเคราะห์การถดถอยในการพยากรณ์ตัวแปรตาม y เมื่อตัวแปรอิสระมีทั้งตัวแปรเชิงคุณภาพและเชิงปริมาณ และมีอันตรกิริยาระหว่างตัวแปรอิสระ พร้อมทั้งแสดงการทดสอบสมมติฐาน และอธิบายความหมายของสมการที่ได้ได้
Task 1: Salary with GPA and Gender
Set working directory (Folder ที่มีไฟล์ข้อมูล)ด้วยคำสั่ง setwd()
และ import data ด้วยคำสั่ง read.csv()
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
dt <- read.csv("SalaryAndGPA.csv")
dt## Obs Salary GPA Gender
## 1 1 49.4 3.0 M
## 2 2 39.0 2.5 M
## 3 3 41.6 1.0 F
## 4 4 54.6 3.5 M
## 5 5 72.8 3.5 F
## 6 6 65.0 3.0 F
## 7 7 54.6 2.0 F
## 8 8 28.6 1.5 M
## 9 9 33.8 2.0 M
## 10 10 44.2 1.5 F#Check data structure
str(dt)## 'data.frame': 10 obs. of 4 variables:
## $ Obs : int 1 2 3 4 5 6 7 8 9 10
## $ Salary: num 49.4 39 41.6 54.6 72.8 65 54.6 28.6 33.8 44.2
## $ GPA : num 3 2.5 1 3.5 3.5 3 2 1.5 2 1.5
## $ Gender: chr "M" "M" "F" "M" ...เปลี่ยนลักษณะของตัวแปร Gender จาก chrเป็น
factor ด้วยคำสั่ง as.factor()
dt$Gender <- as.factor(dt$Gender)
str(dt)## 'data.frame': 10 obs. of 4 variables:
## $ Obs : int 1 2 3 4 5 6 7 8 9 10
## $ Salary: num 49.4 39 41.6 54.6 72.8 65 54.6 28.6 33.8 44.2
## $ GPA : num 3 2.5 1 3.5 3.5 3 2 1.5 2 1.5
## $ Gender: Factor w/ 2 levels "F","M": 2 2 1 2 1 1 1 2 2 1Explore ข้อมูลด้วยสถิติเชิงพรรณนา โดนการสร้างแผนภาพ
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.1#Create scatter plot of Salary and GPA by Gender
ggplot(dt,aes(x=GPA,y=Salary,group=Gender))+
geom_point(aes(col=Gender))
#Create boxplot of Salary by Gender
ggplot(dt,aes(x=Gender,y=Salary))+
geom_boxplot()
#Summary Salary by Gender
aggregate(Salary~Gender,data=dt,FUN=summary)## Gender Salary.Min. Salary.1st Qu. Salary.Median Salary.Mean Salary.3rd Qu.
## 1 F 41.60 44.20 54.60 55.64 65.00
## 2 M 28.60 33.80 39.00 41.08 49.40
## Salary.Max.
## 1 72.80
## 2 54.60กำหนดให้ \(Y\) แทน ตัวแปรตาม Salary ( in thousand dollars)\ และ \(X_{1}\) เป็นตัวแปรอิสระแทน Grade point average (GPA) และ \(X_{2}\)เป็นตัวแปรหุ่นแทน Gender ดังนี้ \[\begin{equation} X_{2}= \begin{cases} 1, & \text{Female} \\ 0, & \text{Male} \end{cases} \end{equation}\]
y <- dt$Salary
x1 <- dt$GPA
x2 <- c() #create null vector for x2
n <- length(y) #number of observations
for(i in 1:n){
if(dt$Gender[i] == "F") x2[i]=1
else x2[i] =0
}
x2## [1] 0 0 1 0 1 1 1 0 0 1ทำการวิเคราะห์ข้อมูลข้างต้นโดยใช้การวิเคราะห์การถดถอยที่มีตัวแบบ ดังนี้ \[\begin{equation} Y = \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{1}X_{2}+\epsilon \end{equation}\]
#Fit linear regression
fit1 <- lm(y~x1+x2+x1*x2)
summary(fit1)##
## Call:
## lm(formula = y ~ x1 + x2 + x1 * x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5093 -0.7649 0.2872 1.2123 1.5600
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.2800 2.9510 2.467 0.04866 *
## x1 13.5200 1.1359 11.903 2.13e-05 ***
## x2 20.2921 3.6034 5.631 0.00134 **
## x1:x2 -0.7619 1.4284 -0.533 0.61294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.796 on 6 degrees of freedom
## Multiple R-squared: 0.9887, Adjusted R-squared: 0.983
## F-statistic: 174.3 on 3 and 6 DF, p-value: 3.178e-06anova(fit1)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 860.29 860.29 266.7213 3.356e-06 ***
## x2 1 825.67 825.67 255.9885 3.786e-06 ***
## x1:x2 1 0.92 0.92 0.2845 0.6129
## Residuals 6 19.35 3.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ทำการวิเคราะห์ข้อมูลข้างต้นโดยใช้การวิเคราะห์การถดถอยที่มีตัวแบบดังนี้
\[\begin{equation}
Y = \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\epsilon
\end{equation}\]
fit2 <- lm(y~x1+x2)
summary(fit2)
anova(fit2)Task 2: ปัจจัยที่ส่งผลต่อคะแนนสอบวิชาคณิตศาสตร์
ทำการศึกษาตัวแปรที่คาดว่าจะส่งผลต่อคะแนนสอบวิชาคณิตศาสตร์ของนักเรียนโรงเรียนแห่งหนึ่ง โดยพิจารณาตัวแปรที่เกี่ยวข้องคือ วิธีการเรียนการสอนที่แตกต่างกัน โดยแบ่งเป็น 4 กลุ่มคือ วิธี A, B, C และ D คะแนนPretest และคะแนนความฉลาดทางอารมณ์(EQ) ทำการสุ่มเก็บข้อมูลจาก 36 หน่วยตัวอย่าง
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
dt2 <- read.csv("ScoreAndEQ.csv")
head(dt2) #show the first six rows## Obs Y Method Pretest EQ
## 1 1 61 B 25 75
## 2 2 78 D 26 79
## 3 3 47 A 20 70
## 4 4 65 B 23 74
## 5 5 45 B 21 69
## 6 6 106 D 40 80#Check data structure
str(dt2)## 'data.frame': 36 obs. of 5 variables:
## $ Obs : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Y : int 61 78 47 65 45 106 120 49 45 62 ...
## $ Method : chr "B" "D" "A" "B" ...
## $ Pretest: int 25 26 20 23 21 40 45 24 26 25 ...
## $ EQ : int 75 79 70 74 69 80 85 82 72 71 ...กำหนดให้ \(Y\) แทน คะแนนสอบวิชาคณิตศาสตร์ (เต็ม 150 คะแนน) และกำหนดตัวแปรหุ่น \(X_{1}, X_{2}, X_{3}\) แทนตัวแปรอิสระวิธีการเรียนการสอนที่มี 4 กลุ่ม ดังนี้
\[\begin{equation} X_{1}= \begin{cases} 1, & \text{วิธีการสอนแบบ A} \\ 0, & \text{อื่นๆ} \end{cases} , X_{2} = \begin{cases} 1, & \text{วิธีการสอนแบบ B} \\ 0, & \text{อื่นๆ} \end{cases} ,X_{3}= \begin{cases} 1, & \text{วิธีการสอนแบบ C} \\ 0, & \text{อื่นๆ} \end{cases} \end{equation}\] โดยให้วิธีการสอนแบบ D เป็นกลุ่มอ้างอิง และให้ \(X_{4}\) แทนคะแนน Pretest และให้ \(X_{5}\) เป็นตัวแปรอิสระแทนคะแนนความฉลาดทางอารมณ์ (E.Q.)
# TO DO ทำการวิเคราะห์ข้อมูลข้างต้นด้วยสถิติเชิงพรรณา (การนำเสนอด้วยตัวเลข และการนำเสนอด้วยแผนภาพ)
# TO DO ทำการวิเคราะห์ข้อมูลข้างต้นโดยใช้การวิเคราะห์การถดถอยที่มีตัวแบบ ดังนี้ \[\begin{equation} Y = \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{3}X_{3}+\beta_{4}X_{4}+\beta_{5}X_{5}+\epsilon \end{equation}\]
# TO DO ปฏิบัติการ: การคัดเลือกตัวแปรอิสระ
นักศึกษาสามารถแสดงการคัดเลือกตัวแปรในการวิเคราะห์การถดถอยที่มีตัวแปรอิสระหลายตัว โดยใช้วิธีการ Forward selection, Backward elimination, และ Stepwise regression ในโปรแกรม R ได้
Task 1: VO2 Oxygen consumption
Import ข้อมูล V02Data.csvและเปลี่ยนชื่อตัวแปร
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
dtv <- read.csv("VO2Data.csv")
head(dtv) #show the first six rows## ID V02 Weight Age Gender HeartRate BodyTemp Height
## 1 1 2.1799 177.0 28 Female 164 37.10 181.17
## 2 2 3.7206 204.1 22 Female 163 37.17 181.12
## 3 3 3.3891 161.0 27 Female 170 37.14 172.65
## 4 4 2.3899 159.3 39 Female 150 37.40 169.91
## 5 5 2.2641 144.5 23 Female 165 37.03 163.21
## 6 6 3.7882 198.3 37 Female 154 37.44 187.66#Check data structure
str(dtv)## 'data.frame': 23 obs. of 8 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ V02 : num 2.18 3.72 3.39 2.39 2.26 ...
## $ Weight : num 177 204 161 159 144 ...
## $ Age : int 28 22 27 39 23 37 25 20 28 39 ...
## $ Gender : chr "Female" "Female" "Female" "Female" ...
## $ HeartRate: int 164 163 170 150 165 154 160 166 150 170 ...
## $ BodyTemp : num 37.1 37.2 37.1 37.4 37 ...
## $ Height : num 181 181 173 170 163 ...y = dtv$V02
x1 = dtv$Weight
x2 = dtv$Age
x4= dtv$HeartRate
x5 = dtv$BodyTemp
x6 = dtv$Heightกำหนดตัวแปรหุ่นใหเกับตัวแปร Gender โดยกำหนดให้เป็นตัวแปร x3 มีค่า = 1 หากเป็น Male และ มีค่า = 0 ถ้าเป็น Female
# dummy coding
x3= c()
for(i in 1:length(dtv$Gender)){
if(dtv$Gender[i]=="Male") x3[i] = 1
else x3[i] = 0
}
x3## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1สร้างตัวแบบการถดถอยที่มีเฉพาะค่าคงที่ และตัวแบบการถดถอยที่มีตัวแปรอิสระ x1 x2 x3 และ x4
# create models
null=lm(y~1, data=dtv)
full=lm(y~x1+x2+x3+x4, data=dtv)ตรวจสอบ multicollinearity
#need to install package "car"
library(car)## Loading required package: carDatavif(full)## x1 x2 x3 x4
## 1.184108 1.055938 1.071945 1.175913#cross-check vif values; vif=1/(1-Rj2)
fitx1 <- lm(x1~x2+x3+x4, data=dtv)
summary(fitx1)##
## Call:
## lm(formula = x1 ~ x2 + x3 + x4, data = dtv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.370 -22.841 -0.861 27.083 40.180
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 437.8522 162.8599 2.689 0.0145 *
## x2 -0.6170 0.9890 -0.624 0.5401
## x3 18.6634 17.7004 1.054 0.3049
## x4 -1.5491 0.9659 -1.604 0.1253
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.41 on 19 degrees of freedom
## Multiple R-squared: 0.1555, Adjusted R-squared: 0.02214
## F-statistic: 1.166 on 3 and 19 DF, p-value: 0.3487fitx2 <- lm(x2~x1+x3+x4, data=dtv)
summary(fitx2)##
## Call:
## lm(formula = x2 ~ x1 + x3 + x4, data = dtv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.0118 -5.2951 -0.8204 5.0059 11.1302
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 69.54565 40.94024 1.699 0.106
## x1 -0.03254 0.05215 -0.624 0.540
## x3 2.06007 4.15501 0.496 0.626
## x4 -0.21321 0.23123 -0.922 0.368
##
## Residual standard error: 6.524 on 19 degrees of freedom
## Multiple R-squared: 0.05297, Adjusted R-squared: -0.09656
## F-statistic: 0.3543 on 3 and 19 DF, p-value: 0.7866แสดงการคัดเลือกตัวแปรอิสระด้วยวิธี Forward selection, Backward elimination และ Stepwise regression
#forward selection
step(null, scope=list(lower=null, upper=full),direction="forward")
#backward elimination
step(full,data=dtv,direction="backward")
# stepwise
step(full, direction="both", data=dtv)ปฏิบัติการ: การตรวจสอบความเหมาะสมของตัวแบบ
นักศึกษาสามารถ
- แสดงการตรวจสอบค่าส่วนเหลือว่าเป็นไปตามข้อสมมติของการวิเคราะห์การถดถอยได้
- แสดงการแก้ไขปัญหาอันเกิดจากค่าส่วนเหลือไม่เป็นไปตามข้อสมมติของการวิเคราะห์การถดถอยได้
Task 1: The Rocket Propellant Data
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
dtr <- read.csv("TheRocketData.csv")
str(dtr)## 'data.frame': 20 obs. of 3 variables:
## $ Observation : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Sher.Strength : num 2159 1678 2316 2061 2208 ...
## $ Age.of.Propellant: num 15.5 23.8 8 17 5.5 ...head(dtr)## Observation Sher.Strength Age.of.Propellant
## 1 1 2158.70 15.50
## 2 2 1678.15 23.75
## 3 3 2316.00 8.00
## 4 4 2061.30 17.00
## 5 5 2207.50 5.50
## 6 6 1708.30 19.001.ทำการวิเคราะห์การถดถอยเชิงเส้นอย่างง่ายด้วยตัวแบบ \[\begin{equation} y = \beta_{0}+\beta_{1}x_{1}+\epsilon \end{equation}\]
x <- dtr$Sher.Strength
y <- dtr$Age.of.Propellant
fit1 <- lm(y ~ x)
summary(fit1)
anova(fit1)2.ทำการตรวจสอบความเหมาะสมของตัวแบบด้วยการวิเคราะห์ค่าส่วนเหลือ ด้วยการตรวจสอบด้วยแผนภาพและคำนวณค่าสถิติทดสอบที่เกี่ยวข้อง พร้อมบันทึกค่าผลลัพธ์ต่างๆที่ได้ลงในตาราง
res <- resid(fit1)
#residuals diagnostic plots
par(mfrow=c(2,2))#create 2x2 plotting layout
plot(fit1)
#Kolmogorov-Smirnov test in R
ks.test(res, "pnorm", mean=mean(res), sd=sd(res))
#Shapiro-Wilk test
shapiro.test(res)
# Breusch-Pagan Test; need to install package "lmtest"
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(fit1)
#Durbin-Watson test; need to install package "lmtest"
library(lmtest)
dwtest(fit1)- ใช้การแปลงด้วยวิธี BoxCox transformation
#BoxCoxTransformation
library(MASS)
par(mfrow=c(1,1))#set plot layout
boxcox(fit1,lambda=seq(-3,3,by=0.05))
4.ทำการวิเคราะห์การถดถอยเชิงเส้นอย่างง่ายด้วยตัวแบบ \[\begin{equation} \sqrt{y} = \beta_{0}+\beta_{1}x_{1}+\epsilon \end{equation}\] และตรวจสอบความเหมาะสมของตัวแบบด้วยการวิเคราะห์ค่าส่วนเหลือ ด้วยการตรวจสอบด้วยแผนภาพและคำนวณค่าสถิติทดสอบที่เกี่ยวข้อง พร้อมบันทึกค่าผลลัพธ์ต่างๆที่ได้ลงในตาราง
# TO DO
fit2 <- lm(sqrt(y) ~ x)Task 2: Olympics Athletes
The data in this task represents the level of cholesterol(in milligrams per liter) and average daily intake of saturated fat (in milligrams) for a sample of 10 Olympics athletes.
setwd("C:\\Users\\ASUS\\Google Drive\\208329\\Term1_2566\\lab-data")
dt <- read.csv("TheAthleteData.csv")
head(dt) #show the first six rows## Athlete Fat.Intake Cholesterol
## 1 1 1290 1182
## 2 2 1350 1172
## 3 3 1470 1264
## 4 4 1600 1493
## 5 5 1710 1571
## 6 6 1840 1711#Check data structure
str(dt)## 'data.frame': 10 obs. of 3 variables:
## $ Athlete : int 1 2 3 4 5 6 7 8 9 10
## $ Fat.Intake : int 1290 1350 1470 1600 1710 1840 1980 2230 2400 2930
## $ Cholesterol: int 1182 1172 1264 1493 1571 1711 1804 1840 1956 1954- Plot a scatter diagram. Does it seem likely that a straight-line model will be adequate?
# TO DO - Fit the straight-line model. Compute the summary statistics and the residual plots. What are your conclusions regarding model adequacy?
# TO DO - Identify an appropriate transformed model for these data. Fit this model to the data and conduct the usual tests of model adequacy.
# TO DO - Write down the summary statistics from the two models in the lab sheet