-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlab4_script.R
91 lines (71 loc) · 2.81 KB
/
lab4_script.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# x - 1st set
# y - 2nd set
# z - 3rd set
mydata = read.csv("Data_21var.csv")
x = mydata[["Spec_Length"]]
y = mydata[["Energy"]]
z = mydata[["Passengers"]]
##### Draw scatter plots
plot(x, y, xlab="x", ylab="y", main="Correlation between X and Y") # scatter plot
# Add trend line to the plot
lm.out = lm(y ~ x)
abline(lm.out, col="red")
# Add 95% confidence intervals for the trend line
newx_xy = seq(min(x),max(x),by = 0.05)
conf_interval = predict(lm.out, newdata=data.frame(x=newx_xy), interval="confidence",
level = 0.95)
#cat("Condidence interval for x-y", conf_interval)
lines(newx_xy, conf_interval[,2], col="blue", lty=2)
lines(newx_xy, conf_interval[,3], col="blue", lty=2)
plot(x, z, xlab="x", ylab="z", main="Correlation between X and Z")
lm.out = lm(z ~ x)
abline(lm.out, col="red")
newx_xz = seq(min(x),max(x),by=0.05)
conf_interval = predict(lm.out, newdata=data.frame(x=newx_xz), interval="confidence",
level = 0.95)
#cat("Condidence interval for x-z", conf_interval)
lines(newx_xz, conf_interval[,2],col="blue", lty=2)
lines(newx_xz, conf_interval[,3],col="blue", lty=2)
plot(y, z, xlab="y", ylab="z", main="Correlation between Y and Z")
lm.out = lm(z ~ y)
abline(lm.out, col="red")
newy_yz = seq(min(y), max(y), by=0.05)
conf_interval = predict(lm.out, newdata=data.frame(y=newy_yz), interval="confidence",
level = 0.95)
#cat("Condidence interval for y-z", conf_interval)
lines(newy_yz, conf_interval[,2],col="blue", lty=2)
lines(newy_yz, conf_interval[,3],col="blue", lty=2)
#plot(z, y, xlab="z", ylab="y", main="Correlation between Z and Y")
#lm.out = lm(y ~ z)
#abline(lm.out, col="red")
#newy_zy = seq(min(z), max(z), by=0.05)
#conf_interval = predict(lm.out, newdata=data.frame(z=newy_zy), interval="confidence",
# level = 0.95)
#cat("Condidence interval for y-z", conf_interval)
#lines(newy_zy, conf_interval[,2],col="blue", lty=2)
#lines(newy_zy, conf_interval[,3],col="blue", lty=2)
##### Calcucate correlation coefficients
cor_xy = cor.test(x, y)
cor_xy
cor_yx = cor.test(y, x)
cor_yx
cor_xz = cor.test(x, z)
cor_xz
cor_yz = cor.test(y, z)
cor_yz
# For non-rejection region calculation use function 'qt(p, df)' to know the
# quantiles of the t distribution. Don't forget to set right p and df parameters
alpha = 0.05
critical_value_l = qt(alpha/2, cor_xy$parameter) # = -2.039513
critical_value_r = qt(1-alpha/2, cor_xy$parameter) # = 2.039513
##### Partial correlation coefficients
source("F:/dev/R/Labworks/pcor.R")
pcor.test(x, y, z)
pcor.test(x, z, y)
#pcor.test(y, x, z)
pcor.test(y, z, x)
# Multiple correlation coefficient
# r=sqrt((rxy^2+rxz^2-2*rxy*rxz*ryz)/(1-ryz^2));
rxyz = sqrt((cor_xy$estimate^2 + cor_xz$estimate^2 - 2*cor_xy$estimate * cor_xz$estimate * cor_yz$estimate) /
(1 - cor_yz$estimate^2))
rxyz