RAnalytix

Data Analysis using R Statistical Package

E-Content For the courses AEC - PGS 504 , Basic Laboratory Techniques , PAJANCOA & STAT 522, Data Analysis using R Statistical Packages , KAU

By

Mr.K.C.Ayyoob, Assistant Professor (Agricultural Statistics)
KAU, Kerala-680656

Dr.S.Saravanan, Associate Professor (Computer Science)
PAJANCOA & RI, Puducherry-609603

R programming language is a free and open-source software used for all kinds of data analysis, statistics, and visualization projects. These notes are an introduction R software for an introductory and advanced level statistics course.

R comes with several built-in data sets, which are generally used as demo data for demonstrating R functions. Majority of the examples here are shown using the built in data sets in R

R - Generating Frequency Tables

Syntax

table(x)

# The R data set iris is used to generate the table

> iris$Sepal.Length

[1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4 5.1

[19] 5.7 5.1 5.4 5.1 4.6 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5 4.9 5.0

[37] 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8 5.1 4.6 5.3 5.0 7.0 6.4 6.9 5.5

[55] 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8 6.2 5.6 5.9 6.1

[73] 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4 6.0 6.7 6.3 5.6 5.5

[91] 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3

[109] 6.7 7.2 6.5 6.4 6.8 5.7 5.8 6.4 6.5 7.7 7.7 6.0 6.9 5.6 7.7 6.3 6.7 7.2

[127] 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7 6.3 6.4 6.0 6.9 6.7 6.9 5.8 6.8

[145] 6.7 6.7 6.3 6.5 6.2 5.9

> table(iris$Sepal.Length)

4.3 4.4 4.5 4.6 4.7 4.8 4.9 5 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6 6.1 6.2

1 3 1 4 2 5 6 10 9 4 1 6 7 6 8 7 3 6 6 4

6.3 6.4 6.5 6.6 6.7 6.8 6.9 7 7.1 7.2 7.3 7.4 7.6 7.7 7.9

9 7 5 2 8 3 4 1 1 3 1 1 1 4 1

R - Summary

> PlantGrowth

weight group

1 4.17 ctrl

2 5.58 ctrl

3 5.18 ctrl

4 6.11 ctrl

5 4.50 ctrl

6 4.61 ctrl

7 5.17 ctrl

8 4.53 ctrl

9 5.33 ctrl

10 5.14 ctrl

11 4.81 trt1

12 4.17 trt1

13 4.41 trt1

14 3.59 trt1

15 5.87 trt1

16 3.83 trt1

17 6.03 trt1

18 4.89 trt1

19 4.32 trt1

20 4.69 trt1

21 6.31 trt2

22 5.12 trt2

23 5.54 trt2

24 5.50 trt2

25 5.37 trt2

26 5.29 trt2

27 4.92 trt2

28 6.15 trt2

29 5.80 trt2

30 5.26 trt2

> summary(PlantGrowth)

weight group

Min. :3.590 ctrl:10

1st Qu.:4.550 trt1:10

Median :5.155 trt2:10

Mean :5.073

3rd Qu.:5.530

Max. :6.310

F Test

To Compare Variances of Two Variables Whether They Are Equal or Not ?

Syntax

var.test(x, y)

> x <- rnorm(50, mean = 0, sd = 2)

> y <- rnorm(30, mean = 1, sd = 1)

> x

[1] -2.38674249 3.02334459 3.00853003 2.80954691 2.27742129 1.00027789

[7] 2.99106354 1.33254057 1.17546519 1.50990704 3.76452844 0.84384436

[13] -0.30968668 -2.82906475 0.44599269 2.23038352 -0.10133583 5.75598167

[19] -1.35609801 -1.73875230 -1.58411426 0.92637318 -0.37979074 -1.49921103

[25] -1.94038175 0.42922356 3.30565031 2.93337499 1.75370754 -3.23368791

[31] 2.53329186 0.44038469 -2.60403946 1.51146349 1.71777220 -0.53925843

[37] 0.51733394 -1.05068354 -2.24446839 0.81881607 1.04579992 -3.56237388

[43] 2.34905136 1.72077751 1.77653050 0.20816790 2.18319760 2.98693591

[49] 0.71411986 -0.01017982

> y

[1] -0.34304280 0.73666863 0.12709703 0.60898492 0.43953553 -0.37665983

[7] 1.17677697 2.17214972 -0.48774027 1.75231874 3.02595490 0.50897167

[13] -0.80756930 2.66499336 2.16658165 -0.27310823 -0.02581134 1.50764241

[19] 3.72589277 3.09258436 0.26652243 0.93384059 1.14415832 2.04060078

[25] 0.72735389 1.49763470 1.91437012 1.07630279 1.88257215 0.35829262

> var.test(x, y)

F test to compare two variances

data: x and y

F = 3.0541, num df = 49, denom df = 29, p-value = 0.001914

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

1.534436 5.745980

sample estimates:

ratio of variances

3.054071

Correlation

> head(Orange)

Tree age circumference

1 1 118 30

2 1 484 58

3 1 664 87

4 1 1004 115

5 1 1231 120

6 1 1372 142

---

> cor(Orange$circumference,Orange$age)

[1] 0.9135189

---

Note:

By Default, it Uses method = "pearson"

For More Details, See Below

---

> longley

GNP.deflator GNP Unemployed Armed.Forces Population Year Employed

1947 83.0 234.289 235.6 159.0 107.608 1947 60.323

1948 88.5 259.426 232.5 145.6 108.632 1948 61.122

1949 88.2 258.054 368.2 161.6 109.773 1949 60.171

1950 89.5 284.599 335.1 165.0 110.929 1950 61.187

1951 96.2 328.975 209.9 309.9 112.075 1951 63.221

1952 98.1 346.999 193.2 359.4 113.270 1952 63.639

1953 99.0 365.385 187.0 354.7 115.094 1953 64.989

1954 100.0 363.112 357.8 335.0 116.219 1954 63.761

1955 101.2 397.469 290.4 304.8 117.388 1955 66.019

1956 104.6 419.180 282.2 285.7 118.734 1956 67.857

1957 108.4 442.769 293.6 279.8 120.445 1957 68.169

1958 110.8 444.546 468.1 263.7 121.950 1958 66.513

1959 112.6 482.704 381.3 255.2 123.366 1959 68.655

1960 114.2 502.601 393.1 251.4 125.368 1960 69.564

1961 115.7 518.173 480.6 257.2 127.852 1961 69.331

1962 116.9 554.894 400.7 282.7 130.081 1962 70.551

---

> options("width"=120)

---

> cor(longley, method = "spearman")

GNP.deflator GNP Unemployed Armed.Forces Population Year Employed

GNP.deflator 1.0000000 0.9970588 0.6647059 0.2205882 0.9970588 0.9970588 0.9823529

GNP 0.9970588 1.0000000 0.6382353 0.2235294 0.9941176 0.9941176 0.9852941

Unemployed 0.6647059 0.6382353 1.0000000 -0.3411765 0.6852941 0.6852941 0.5647059

Armed.Forces 0.2205882 0.2235294 -0.3411765 1.0000000 0.2264706 0.2264706 0.2264706

Population 0.9970588 0.9941176 0.6852941 0.2264706 1.0000000 1.0000000 0.9764706

Year 0.9970588 0.9941176 0.6852941 0.2264706 1.0000000 1.0000000 0.9764706

Employed 0.9823529 0.9852941 0.5647059 0.2264706 0.9764706 0.9764706 1.0000000

---

> cor(longley, method = "kendall")

GNP.deflator GNP Unemployed Armed.Forces Population Year Employed

GNP.deflator 1.00000000 0.9833333 0.4500000 0.03333333 0.9833333 0.9833333 0.9166667

GNP 0.98333333 1.0000000 0.4333333 0.05000000 0.9666667 0.9666667 0.9333333

Unemployed 0.45000000 0.4333333 1.0000000 -0.21666667 0.4666667 0.4666667 0.3666667

Armed.Forces 0.03333333 0.0500000 -0.2166667 1.00000000 0.0500000 0.0500000 0.0500000

Population 0.98333333 0.9666667 0.4666667 0.05000000 1.0000000 1.0000000 0.9000000

Year 0.98333333 0.9666667 0.4666667 0.05000000 1.0000000 1.0000000 0.9000000

Employed 0.91666667 0.9333333 0.3666667 0.05000000 0.9000000 0.9000000 1.0000000

---

> cor(longley, method = "pearson")

GNP.deflator GNP Unemployed Armed.Forces Population Year Employed

GNP.deflator 1.0000000 0.9915892 0.6206334 0.4647442 0.9791634 0.9911492 0.9708985

GNP 0.9915892 1.0000000 0.6042609 0.4464368 0.9910901 0.9952735 0.9835516

Unemployed 0.6206334 0.6042609 1.0000000 -0.1774206 0.6865515 0.6682566 0.5024981

Armed.Forces 0.4647442 0.4464368 -0.1774206 1.0000000 0.3644163 0.4172451 0.4573074

Population 0.9791634 0.9910901 0.6865515 0.3644163 1.0000000 0.9939528 0.9603906

Year 0.9911492 0.9952735 0.6682566 0.4172451 0.9939528 1.0000000 0.9713295

Employed 0.9708985 0.9835516 0.5024981 0.4573074 0.9603906 0.9713295 1.0000000

---

For More Details:

> ?cor

ANOVA Without Interactions Between Factors (N, P & K)

Syntax

aov(<dependent variable>~<independent variable1>+<independent variable2>+..., data=<data frame>)

---

> summary(aov(yield ~ block + N+P+K, data=npk))

Df Sum Sq Mean Sq F value Pr(>F)

block 5 343.3 68.66 4.288 0.01272 *

N 1 189.3 189.28 11.821 0.00366 **

P 1 8.4 8.40 0.525 0.47999

K 1 95.2 95.20 5.946 0.02767 *

Residuals 15 240.2 16.01

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANOVA With Interactions Between Factors (N, P & K)

Syntax

aov(<dependent variable>~<independent variable1>*<independent variable2>*..., data=<data frame>)

---

> summary(aov(yield ~ block + N*P*K, data=npk))

Df Sum Sq Mean Sq F value Pr(>F)

block 5 343.3 68.66 4.447 0.01594 *

N 1 189.3 189.28 12.259 0.00437 **

P 1 8.4 8.40 0.544 0.47490

K 1 95.2 95.20 6.166 0.02880 *

N:P 1 21.3 21.28 1.378 0.26317

N:K 1 33.1 33.13 2.146 0.16865

P:K 1 0.5 0.48 0.031 0.86275

Residuals 12 185.3 15.44

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANCOVA (Analyisis of Covariance)

> head(mtcars)

mpg cyl disp hp drat wt qsec vs am gear carb

Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4

Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4

Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1

Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1

Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2

Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1

> summary(aov(mpg~hp+am,data=mtcars))

Df Sum Sq Mean Sq F value Pr(>F)

hp 1 678.4 678.4 80.15 7.63e-10 ***

am 1 202.2 202.2 23.89 3.46e-05 ***

Residuals 29 245.4 8.5

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Simple Linear Regression (lm - Linear Model)

Syntax

summary(lm(<dependent variable>~<independent variable>, data=<data frame>))

---

> summary(lm(circumference~age, data=Orange))

Call:

lm(formula = circumference ~ age, data = Orange)

Residuals:

Min 1Q Median 3Q Max

-46.310 -14.946 -0.076 19.697 45.111

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 17.399650 8.622660 2.018 0.0518 .

age 0.106770 0.008277 12.900 1.93e-14 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.74 on 33 degrees of freedom

Multiple R-squared: 0.8345, Adjusted R-squared: 0.8295

F-statistic: 166.4 on 1 and 33 DF, p-value: 1.931e-14

Multiple Linear Regression (lm - Linear Model)

Syntax

summary(lm(<dependent variable>~<independent variable1 + independent variable2 + ... >, data=<data frame>))

---

> summary(lm(VanKilled~PetrolPrice + kms, data=Seatbelts))

Call:

lm(formula = VanKilled ~ PetrolPrice + kms, data = Seatbelts)

Residuals:

Min 1Q Median 3Q Max

-9.0120 -2.3903 0.3601 2.3563 6.5804

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 2.102e+01 2.003e+00 10.495 < 2e-16 ***

PetrolPrice -3.411e+01 2.025e+01 -1.684 0.0938 .

kms -5.622e-04 8.393e-05 -6.699 2.34e-10 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.147 on 189 degrees of freedom

Multiple R-squared: 0.2592, Adjusted R-squared: 0.2513

F-statistic: 33.06 on 2 and 189 DF, p-value: 4.888e-13

Maximum Likelihood Estimation (MLE)

> library("stats4")

> x <- rnorm(100, mean = 3, sd = 2)

> functiontocal <- function(mu, sigma)

{

R = dnorm(x, mu, sigma)

-sum(log(R))

}

> summary(mle(functiontocal, start = list(mu = 1, sigma=1)))

Maximum likelihood estimation

Call:

mle(minuslogl = functiontocal, start = list(mu = 1, sigma = 1))

Coefficients:

Estimate Std. Error

mu 3.243171 0.1797415

sigma 1.797415 0.1270963

-2 log L: 401.0576

Warning messages:

1: In dnorm(x, mu, sigma) : NaNs produced

2: In dnorm(x, mu, sigma) : NaNs produced

...

Econometrics - Logit Model (Logistic Regression / Logit Regression), Probit Model

Download Data From Below Link & Have it in $HOME Folder

https://stats.idre.ucla.edu/stat/data/binary.csv

> mydata=read.csv("binary.csv")

> mydata

admit gre gpa rank

1 0 380 3.61 3

2 1 660 3.67 3

3 1 800 4.00 1

4 1 640 3.19 4

5 0 520 2.93 4

6 1 760 3.00 2

7 1 560 2.98 1

8 0 400 3.08 2

9 1 540 3.39 3

10 0 700 3.92 2

11 0 800 4.00 4

12 0 440 3.22 1

13 1 760 4.00 1

14 0 700 3.08 2

15 1 700 4.00 1

16 0 480 3.44 3

17 0 780 3.87 4

18 0 360 2.56 3

19 0 800 3.75 2

20 1 540 3.81 1

21 0 500 3.17 3

22 1 660 3.63 2

23 0 600 2.82 4

24 0 680 3.19 4

25 1 760 3.35 2

26 1 800 3.66 1

27 1 620 3.61 1

28 1 520 3.74 4

29 1 780 3.22 2

30 0 520 3.29 1

31 0 540 3.78 4

32 0 760 3.35 3

33 0 600 3.40 3

34 1 800 4.00 3

35 0 360 3.14 1

36 0 400 3.05 2

37 0 580 3.25 1

38 0 520 2.90 3

39 1 500 3.13 2

40 1 520 2.68 3

41 0 560 2.42 2

42 1 580 3.32 2

43 1 600 3.15 2

44 0 500 3.31 3

45 0 700 2.94 2

46 1 460 3.45 3

47 1 580 3.46 2

48 0 500 2.97 4

49 0 440 2.48 4

50 0 400 3.35 3

51 0 640 3.86 3

52 0 440 3.13 4

53 0 740 3.37 4

54 1 680 3.27 2

55 0 660 3.34 3

56 1 740 4.00 3

57 0 560 3.19 3

58 0 380 2.94 3

59 0 400 3.65 2

60 0 600 2.82 4

61 1 620 3.18 2

62 0 560 3.32 4

63 0 640 3.67 3

64 1 680 3.85 3

65 0 580 4.00 3

66 0 600 3.59 2

67 0 740 3.62 4

68 0 620 3.30 1

69 0 580 3.69 1

70 0 800 3.73 1

71 0 640 4.00 3

72 0 300 2.92 4

73 0 480 3.39 4

74 0 580 4.00 2

75 0 720 3.45 4

76 0 720 4.00 3

77 0 560 3.36 3

78 1 800 4.00 3

79 0 540 3.12 1

80 1 620 4.00 1

81 0 700 2.90 4

82 0 620 3.07 2

83 0 500 2.71 2

84 0 380 2.91 4

85 1 500 3.60 3

86 0 520 2.98 2

87 0 600 3.32 2

88 0 600 3.48 2

89 0 700 3.28 1

90 1 660 4.00 2

91 0 700 3.83 2

92 1 720 3.64 1

93 0 800 3.90 2

94 0 580 2.93 2

95 1 660 3.44 2

96 0 660 3.33 2

97 0 640 3.52 4

98 0 480 3.57 2

99 0 700 2.88 2

100 0 400 3.31 3

101 0 340 3.15 3

102 0 580 3.57 3

103 0 380 3.33 4

104 0 540 3.94 3

105 1 660 3.95 2

106 1 740 2.97 2

107 1 700 3.56 1

108 0 480 3.13 2

109 0 400 2.93 3

110 0 480 3.45 2

111 0 680 3.08 4

112 0 420 3.41 4

113 0 360 3.00 3

114 0 600 3.22 1

115 0 720 3.84 3

116 0 620 3.99 3

117 1 440 3.45 2

118 0 700 3.72 2

119 1 800 3.70 1

120 0 340 2.92 3

121 1 520 3.74 2

122 1 480 2.67 2

123 0 520 2.85 3

124 0 500 2.98 3

125 0 720 3.88 3

126 0 540 3.38 4

127 1 600 3.54 1

128 0 740 3.74 4

129 0 540 3.19 2

130 0 460 3.15 4

131 1 620 3.17 2

132 0 640 2.79 2

133 0 580 3.40 2

134 0 500 3.08 3

135 0 560 2.95 2

136 0 500 3.57 3

137 0 560 3.33 4

138 0 700 4.00 3

139 0 620 3.40 2

140 1 600 3.58 1

141 0 640 3.93 2

142 1 700 3.52 4

143 0 620 3.94 4

144 0 580 3.40 3

145 0 580 3.40 4

146 0 380 3.43 3

147 0 480 3.40 2

148 0 560 2.71 3

149 1 480 2.91 1

150 0 740 3.31 1

151 1 800 3.74 1

152 0 400 3.38 2

153 1 640 3.94 2

154 0 580 3.46 3

155 0 620 3.69 3

156 1 580 2.86 4

157 0 560 2.52 2

158 1 480 3.58 1

159 0 660 3.49 2

160 0 700 3.82 3

161 0 600 3.13 2

162 0 640 3.50 2

163 1 700 3.56 2

164 0 520 2.73 2

165 0 580 3.30 2

166 0 700 4.00 1

167 0 440 3.24 4

168 0 720 3.77 3

169 0 500 4.00 3

170 0 600 3.62 3

171 0 400 3.51 3

172 0 540 2.81 3

173 0 680 3.48 3

174 1 800 3.43 2

175 0 500 3.53 4

176 1 620 3.37 2

177 0 520 2.62 2

178 1 620 3.23 3

179 0 620 3.33 3

180 0 300 3.01 3

181 0 620 3.78 3

182 0 500 3.88 4

183 0 700 4.00 2

184 1 540 3.84 2

185 0 500 2.79 4

186 0 800 3.60 2

187 0 560 3.61 3

188 0 580 2.88 2

189 0 560 3.07 2

190 0 500 3.35 2

191 1 640 2.94 2

192 0 800 3.54 3

193 0 640 3.76 3

194 0 380 3.59 4

195 1 600 3.47 2

196 0 560 3.59 2

197 0 660 3.07 3

198 1 400 3.23 4

199 0 600 3.63 3

200 0 580 3.77 4

201 0 800 3.31 3

202 1 580 3.20 2

203 1 700 4.00 1

204 0 420 3.92 4

205 1 600 3.89 1

206 1 780 3.80 3

207 0 740 3.54 1

208 1 640 3.63 1

209 0 540 3.16 3

210 0 580 3.50 2

211 0 740 3.34 4

212 0 580 3.02 2

213 0 460 2.87 2

214 0 640 3.38 3

215 1 600 3.56 2

216 1 660 2.91 3

217 0 340 2.90 1

218 1 460 3.64 1

219 0 460 2.98 1

220 1 560 3.59 2

221 0 540 3.28 3

222 0 680 3.99 3

223 1 480 3.02 1

224 0 800 3.47 3

225 0 800 2.90 2

226 1 720 3.50 3

227 0 620 3.58 2

228 0 540 3.02 4

229 0 480 3.43 2

230 1 720 3.42 2

231 0 580 3.29 4

232 0 600 3.28 3

233 0 380 3.38 2

234 0 420 2.67 3

235 1 800 3.53 1

236 0 620 3.05 2

237 1 660 3.49 2

238 0 480 4.00 2

239 0 500 2.86 4

240 0 700 3.45 3

241 0 440 2.76 2

242 1 520 3.81 1

243 1 680 2.96 3

244 0 620 3.22 2

245 0 540 3.04 1

246 0 800 3.91 3

247 0 680 3.34 2

248 0 440 3.17 2

249 0 680 3.64 3

250 0 640 3.73 3

251 0 660 3.31 4

252 0 620 3.21 4

253 1 520 4.00 2

254 1 540 3.55 4

255 1 740 3.52 4

256 0 640 3.35 3

257 1 520 3.30 2

258 1 620 3.95 3

259 0 520 3.51 2

260 0 640 3.81 2

261 0 680 3.11 2

262 0 440 3.15 2

263 1 520 3.19 3

264 1 620 3.95 3

265 1 520 3.90 3

266 0 380 3.34 3

267 0 560 3.24 4

268 1 600 3.64 3

269 1 680 3.46 2

270 0 500 2.81 3

271 1 640 3.95 2

272 0 540 3.33 3

273 1 680 3.67 2

274 0 660 3.32 1

275 0 520 3.12 2

276 1 600 2.98 2

277 0 460 3.77 3

278 1 580 3.58 1

279 1 680 3.00 4

280 1 660 3.14 2

281 0 660 3.94 2

282 0 360 3.27 3

283 0 660 3.45 4

284 0 520 3.10 4

285 1 440 3.39 2

286 0 600 3.31 4

287 1 800 3.22 1

288 1 660 3.70 4

289 0 800 3.15 4

290 0 420 2.26 4

291 1 620 3.45 2

292 0 800 2.78 2

293 0 680 3.70 2

294 0 800 3.97 1

295 0 480 2.55 1

296 0 520 3.25 3

297 0 560 3.16 1

298 0 460 3.07 2

299 0 540 3.50 2

300 0 720 3.40 3

301 0 640 3.30 2

302 1 660 3.60 3

303 1 400 3.15 2

304 1 680 3.98 2

305 0 220 2.83 3

306 0 580 3.46 4

307 1 540 3.17 1

308 0 580 3.51 2

309 0 540 3.13 2

310 0 440 2.98 3

311 0 560 4.00 3

312 0 660 3.67 2

313 0 660 3.77 3

314 1 520 3.65 4

315 0 540 3.46 4

316 1 300 2.84 2

317 1 340 3.00 2

318 1 780 3.63 4

319 1 480 3.71 4

320 0 540 3.28 1

321 0 460 3.14 3

322 0 460 3.58 2

323 0 500 3.01 4

324 0 420 2.69 2

325 0 520 2.70 3

326 0 680 3.90 1

327 0 680 3.31 2

328 1 560 3.48 2

329 0 580 3.34 2

330 0 500 2.93 4

331 0 740 4.00 3

332 0 660 3.59 3

333 0 420 2.96 1

334 0 560 3.43 3

335 1 460 3.64 3

336 1 620 3.71 1

337 0 520 3.15 3

338 0 620 3.09 4

339 0 540 3.20 1

340 1 660 3.47 3

341 0 500 3.23 4

342 1 560 2.65 3

343 0 500 3.95 4

344 0 580 3.06 2

345 0 520 3.35 3

346 0 500 3.03 3

347 0 600 3.35 2

348 0 580 3.80 2

349 0 400 3.36 2

350 0 620 2.85 2

351 1 780 4.00 2

352 0 620 3.43 3

353 1 580 3.12 3

354 0 700 3.52 2

355 1 540 3.78 2

356 1 760 2.81 1

357 0 700 3.27 2

358 0 720 3.31 1

359 1 560 3.69 3

360 0 720 3.94 3

361 1 520 4.00 1

362 1 540 3.49 1

363 0 680 3.14 2

364 0 460 3.44 2

365 1 560 3.36 1

366 0 480 2.78 3

367 0 460 2.93 3

368 0 620 3.63 3

369 0 580 4.00 1

370 0 800 3.89 2

371 1 540 3.77 2

372 1 680 3.76 3

373 1 680 2.42 1

374 1 620 3.37 1

375 0 560 3.78 2

376 0 560 3.49 4

377 0 620 3.63 2

378 1 800 4.00 2

379 0 640 3.12 3

380 0 540 2.70 2

381 0 700 3.65 2

382 1 540 3.49 2

383 0 540 3.51 2

384 0 660 4.00 1

385 1 480 2.62 2

386 0 420 3.02 1

387 1 740 3.86 2

388 0 580 3.36 2

389 0 640 3.17 2

390 0 640 3.51 2

391 1 800 3.05 2

392 1 660 3.88 2

393 1 600 3.38 3

394 1 620 3.75 2

395 1 460 3.99 3

396 0 620 4.00 2

397 0 560 3.04 3

398 0 460 2.63 2

399 0 700 3.65 2

400 0 600 3.89 3

---

> summary(glm(admit ~ gre + gpa + rank, data = mydata, family = binomial(link="logit")))

Call:

glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "logit"),

data = mydata)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.5802 -0.8848 -0.6382 1.1575 2.1732

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -3.449548 1.132846 -3.045 0.00233 **

gre 0.002294 0.001092 2.101 0.03564 *

gpa 0.777014 0.327484 2.373 0.01766 *

rank -0.560031 0.127137 -4.405 1.06e-05 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 499.98 on 399 degrees of freedom

Residual deviance: 459.44 on 396 degrees of freedom

AIC: 467.44

Number of Fisher Scoring iterations: 4

---

> summary(glm(admit ~ gre + gpa + rank, data = mydata, family = binomial(link="probit")))

Call:

glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "probit"),

data = mydata)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.5626 -0.8920 -0.6403 1.1631 2.2097

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -2.0915037 0.6718360 -3.113 0.00185 **

gre 0.0013982 0.0006487 2.156 0.03112 *

gpa 0.4643598 0.1950263 2.381 0.01727 *

rank -0.3317117 0.0745524 -4.449 8.61e-06 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 499.98 on 399 degrees of freedom

Residual deviance: 459.48 on 396 degrees of freedom

AIC: 467.48

Number of Fisher Scoring iterations: 4

Stepwise Regression

> head(swiss)

Fertility Agriculture Examination Education Catholic

Courtelary 80.2 17.0 15 12 9.96

Delemont 83.1 45.1 6 9 84.84

Franches-Mnt 92.5 39.7 5 5 93.40

Moutier 85.8 36.5 12 7 33.77

Neuveville 76.9 43.5 17 15 5.16

Porrentruy 76.1 35.3 9 7 90.57

Infant.Mortality

Courtelary 22.2

Delemont 22.2

Franches-Mnt 20.2

Moutier 20.3

Neuveville 20.6

Porrentruy 26.6

Note:

. (dot) in command represents all variables of the dataset

Example

lm(Fertility ~ .,data=swiss)

means

lm(Fertility ~ Fertility + Agriculture + Examination + Education + Catholic + Infant.Mortality, data=swiss)

> summary(step(lm(Fertility ~ .,data=swiss)))

Start: AIC=190.69

Fertility ~ Agriculture + Examination + Education + Catholic +

Infant.Mortality

Df Sum of Sq RSS AIC

- Examination 1 53.03 2158.1 189.86

<none> 2105.0 190.69

- Agriculture 1 307.72 2412.8 195.10

- Infant.Mortality 1 408.75 2513.8 197.03

- Catholic 1 447.71 2552.8 197.75

- Education 1 1162.56 3267.6 209.36

Step: AIC=189.86

Fertility ~ Agriculture + Education + Catholic + Infant.Mortality

Df Sum of Sq RSS AIC

<none> 2158.1 189.86

- Agriculture 1 264.18 2422.2 193.29

- Infant.Mortality 1 409.81 2567.9 196.03

- Catholic 1 956.57 3114.6 205.10

- Education 1 2249.97 4408.0 221.43

Call:

lm(formula = Fertility ~ Agriculture + Education + Catholic +

Infant.Mortality, data = swiss)

Residuals:

Min 1Q Median 3Q Max

-14.6765 -6.0522 0.7514 3.1664 16.1422

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 62.10131 9.60489 6.466 8.49e-08 ***

Agriculture -0.15462 0.06819 -2.267 0.02857 *

Education -0.98026 0.14814 -6.617 5.14e-08 ***

Catholic 0.12467 0.02889 4.315 9.50e-05 ***

Infant.Mortality 1.07844 0.38187 2.824 0.00722 **

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.168 on 42 degrees of freedom

Multiple R-squared: 0.6993, Adjusted R-squared: 0.6707

F-statistic: 24.42 on 4 and 42 DF, p-value: 1.717e-10

Kolmogorov-Smirnov Test

> x = rnorm(50)

> x

[1] 1.27469059 0.32901751 1.08547795 -1.53967037 -0.68308888 0.15584386

[7] 1.80373688 0.26680885 -0.40074455 1.38210328 -0.33640360 0.62504762

[13] 0.53823228 -1.07797092 1.03937083 -0.24990742 0.40530596 -0.19115800

[19] 2.27244935 0.25668865 -0.08309251 -0.62805580 0.46322678 -1.24595469

[25] -1.60742418 1.12389246 -1.24154784 -1.01526393 1.08072601 0.20509969

[31] -0.09811376 -0.42752724 -0.42445595 -0.43294211 0.69396285 -1.10617623

[37] 0.93869653 -0.33207389 0.52751335 1.07121654 -0.35511242 0.13009504

[43] 0.29571333 0.39063429 0.93494733 1.92913308 0.36175055 1.41471870

[49] 0.17871536 0.23004740

> y = runif(30)

> y

[1] 0.81475955 0.54984034 0.43774511 0.02959697 0.39997636 0.21102408

[7] 0.22048581 0.54838088 0.35041053 0.51468835 0.85241835 0.65948428

[13] 0.84374469 0.38413933 0.79276582 0.96877940 0.90075692 0.59430388

[19] 0.57028071 0.76543719 0.49153308 0.42906595 0.76723953 0.64221234

[25] 0.20141944 0.50051108 0.98139722 0.14993413 0.69127872 0.69971360

> ks.test(x,y)

Two-sample Kolmogorov-Smirnov test

data: x and y

D = 0.41333, p-value = 0.002282

alternative hypothesis: two-sided

Wilcoxon Signed Rank Test

> x = c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30)

> y = c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)

> wilcox.test(x, y, paired=TRUE)

Wilcoxon signed rank test

data: x and y

V = 40, p-value = 0.03906

alternative hypothesis: true location shift is not equal to 0

Mann-Whitney U Test (Two Sample Unpaired Wilcoxon Test)

> x = c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)

> y = c(1.15, 0.88, 0.90, 0.74, 1.21)

> wilcox.test(x, y)

Wilcoxon rank sum test

data: x and y

W = 35, p-value = 0.2544

alternative hypothesis: true location shift is not equal to 0

Kruskal-Wallis Test

> kruskal.test(Ozone ~ Month, data = airquality)

Kruskal-Wallis rank sum test

data: Ozone by Month

Kruskal-Wallis chi-squared = 29.267, df = 4, p-value = 6.901e-06

McNemar's Test

> Performance = matrix(c(794, 86, 150, 570), nrow = 2)

> Performance

[,1] [,2]

[1,] 794 150

[2,] 86 570

> mcnemar.test(Performance)

McNemar's Chi-squared test with continuity correction

data: Performance

McNemar's chi-squared = 16.818, df = 1, p-value = 4.115e-05

Factor Analysis

m1 <- cbind(

v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6),

v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5),

v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6),

v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4),

v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5),

v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4))

> m1

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 1 1 3 3 1 1

[2,] 1 2 3 3 1 1

[3,] 1 1 3 4 1 1

[4,] 1 1 3 3 1 2

[5,] 1 1 3 3 1 1

[6,] 1 1 1 1 3 3

[7,] 1 2 1 1 3 3

[8,] 1 1 1 2 3 3

[9,] 1 2 1 1 3 4

[10,] 1 1 1 1 3 3

[11,] 3 3 1 1 1 1

[12,] 3 4 1 1 1 1

[13,] 3 3 1 2 1 1

[14,] 3 3 1 1 1 2

[15,] 3 3 1 1 1 1

[16,] 4 4 5 5 6 6

[17,] 5 6 4 6 4 5

[18,] 6 5 6 4 5 4

--- 6 Variables Reducing it to 3 Factors

> factanal(m1,3)

Call:

factanal(x = m1, factors = 3)

Uniquenesses:

[1] 0.005 0.101 0.005 0.224 0.084 0.005

Loadings:

Factor1 Factor2 Factor3

[1,] 0.944 0.182 0.267

[2,] 0.905 0.235 0.159

[3,] 0.236 0.210 0.946

[4,] 0.180 0.242 0.828

[5,] 0.242 0.881 0.286

[6,] 0.193 0.959 0.196

Factor1 Factor2 Factor3

SS loadings 1.893 1.886 1.797

Proportion Var 0.316 0.314 0.300

Cumulative Var 0.316 0.630 0.929

The degrees of freedom for the model is 0 and the fit was 0.4755

Principal Component Analysis

> head(USArrests)

Murder Assault UrbanPop Rape

Alabama 13.2 236 58 21.2

Alaska 10.0 263 48 44.5

Arizona 8.1 294 80 31.0

Arkansas 8.8 190 50 19.5

California 9.0 276 91 40.6

Colorado 7.9 204 78 38.7

---

> prcomp(USArrests)

Standard deviations:

[1] 83.732400 14.212402 6.489426 2.482790

Rotation:

PC1 PC2 PC3 PC4

Murder 0.04170432 -0.04482166 0.07989066 -0.99492173

Assault 0.99522128 -0.05876003 -0.06756974 0.03893830

UrbanPop 0.04633575 0.97685748 -0.20054629 -0.05816914

Rape 0.07515550 0.20071807 0.97408059 0.07232502

---

> summary(prcomp(USArrests))

Importance of components:

PC1 PC2 PC3 PC4

Standard deviation 83.7324 14.21240 6.4894 2.48279

Proportion of Variance 0.9655 0.02782 0.0058 0.00085

Cumulative Proportion 0.9655 0.99335 0.9991 1.00000

R - Example For Creating / Testing Time Series Objects

Syntax

ts(data, frequency, start)

start: time of the first observation

end: time of the last observation

data: a vector or matrix of the observed time-series values

frequency: the number of observations per unit of time

Example 1:

x = ts(c(3,4,5,6,8,11,12,13,14), frequency = 4, start = c(2016, 2))

x

Qtr1 Qtr2 Qtr3 Qtr4

2016 3 4 5

2017 6 8 11 12

2018 13 14

is.ts(x)

[1] TRUE

---

Example 2:

x = ts(c(3,4,5,6,8,11,12,13,14), frequency = 12, start = c(2016, 2))

x

Feb Mar Apr May Jun Jul Aug Sep Oct

2016 3 4 5 6 8 11 12 13 14

is.ts(x)

[1] TRUE

---

Example 3:

x = ts(c(3,4,5,6,8,11,12,13,14), frequency = 3, start = c(2016, 2))

x

Time Series:

Start = c(2016, 2)

End = c(2019, 1)

Frequency = 3

[1] 3 4 5 6 8 11 12 13 14

is.ts(x)

[1] TRUE

R - Example For AutoRegressive (AR) Modelling of Time Series & Predicting

sunspot.year

Time Series:

Start = 1700

End = 1988

Frequency = 1

[1] 5.0 11.0 16.0 23.0 36.0 58.0 29.0 20.0 10.0 8.0 3.0 0.0

[13] 0.0 2.0 11.0 27.0 47.0 63.0 60.0 39.0 28.0 26.0 22.0 11.0

[25] 21.0 40.0 78.0 122.0 103.0 73.0 47.0 35.0 11.0 5.0 16.0 34.0

[37] 70.0 81.0 111.0 101.0 73.0 40.0 20.0 16.0 5.0 11.0 22.0 40.0

[49] 60.0 80.9 83.4 47.7 47.8 30.7 12.2 9.6 10.2 32.4 47.6 54.0

[61] 62.9 85.9 61.2 45.1 36.4 20.9 11.4 37.8 69.8 106.1 100.8 81.6

[73] 66.5 34.8 30.6 7.0 19.8 92.5 154.4 125.9 84.8 68.1 38.5 22.8

[85] 10.2 24.1 82.9 132.0 130.9 118.1 89.9 66.6 60.0 46.9 41.0 21.3

[97] 16.0 6.4 4.1 6.8 14.5 34.0 45.0 43.1 47.5 42.2 28.1 10.1

[109] 8.1 2.5 0.0 1.4 5.0 12.2 13.9 35.4 45.8 41.1 30.1 23.9

[121] 15.6 6.6 4.0 1.8 8.5 16.6 36.3 49.6 64.2 67.0 70.9 47.8

[133] 27.5 8.5 13.2 56.9 121.5 138.3 103.2 85.7 64.6 36.7 24.2 10.7

[145] 15.0 40.1 61.5 98.5 124.7 96.3 66.6 64.5 54.1 39.0 20.6 6.7

[157] 4.3 22.7 54.8 93.8 95.8 77.2 59.1 44.0 47.0 30.5 16.3 7.3

[169] 37.6 74.0 139.0 111.2 101.6 66.2 44.7 17.0 11.3 12.4 3.4 6.0

[181] 32.3 54.3 59.7 63.7 63.5 52.2 25.4 13.1 6.8 6.3 7.1 35.6

[193] 73.0 85.1 78.0 64.0 41.8 26.2 26.7 12.1 9.5 2.7 5.0 24.4

[205] 42.0 63.5 53.8 62.0 48.5 43.9 18.6 5.7 3.6 1.4 9.6 47.4

[217] 57.1 103.9 80.6 63.6 37.6 26.1 14.2 5.8 16.7 44.3 63.9 69.0

[229] 77.8 64.9 35.7 21.2 11.1 5.7 8.7 36.1 79.7 114.4 109.6 88.8

[241] 67.8 47.5 30.6 16.3 9.6 33.2 92.6 151.6 136.3 134.7 83.9 69.4

[253] 31.5 13.9 4.4 38.0 141.7 190.2 184.8 159.0 112.3 53.9 37.5 27.9

[265] 10.2 15.1 47.0 93.8 105.9 105.5 104.5 66.6 68.9 38.0 34.5 15.5

[277] 12.6 27.5 92.5 155.4 154.7 140.5 115.9 66.6 45.9 17.9 13.4 29.2

[289] 100.2

arsunspot = ar(sunspot.year)

arsunspot

Call:

ar(x = sunspot.year)

Coefficients:

1 2 3 4 5 6 7 8

1.1305 -0.3524 -0.1745 0.1403 -0.1358 0.0963 -0.0556 0.0076

9

0.1941

Order selected 9 sigma^2 estimated as 267.5

predict(arsunspot, n.ahead = 25)

$pred

Time Series:

Start = 1989

End = 2013

Frequency = 1

[1] 135.25933 148.09051 133.98476 106.61344 71.21921 40.84057 18.70100

[8] 11.52416 27.24208 56.99888 87.86705 107.62926 111.05437 98.05484

[15] 74.84085 48.80128 27.65441 18.15075 23.15355 40.04723 61.95906

[22] 80.79092 90.11420 87.44131 74.42284

$se

Time Series:

Start = 1989

End = 2013

Frequency = 1

[1] 16.35519 24.68467 28.95653 29.97401 30.07714 30.15629 30.35971 30.58793

[9] 30.71100 30.74276 31.42565 32.96467 34.48910 35.33601 35.51890 35.52034

[17] 35.65505 35.90628 36.07084 36.08139 36.16818 36.56324 37.16527 37.64820

[25] 37.83954

Note: $pred - Prediction, $se - Standard Error

R - Example For MA (Moving Average)

Syntax

SMA(x, n)

x: Price, Volume, etc... Data Series

n: Number of Periods To Average Over

Note: SMA - Simple Moving Average Calculates the Average of the Series Over the Past n Observations

Example

library("TTR")

data("ttrc")

x = ttrc$Close[1:50]

x

[1] 3.08 3.11 3.09 3.10 3.11 3.16 3.22 3.23 3.32 3.30 3.32 3.26 3.16 3.25 3.29

[16] 3.31 3.32 3.37 3.42 3.42 3.44 3.41 3.40 3.43 3.46 3.39 3.43 3.47 3.48 3.50

[31] 3.58 3.54 3.51 3.49 3.48 3.46 3.46 3.45 3.50 3.49 3.51 3.55 3.48 3.49 3.44

[46] 3.40 3.40 3.39 3.42 3.37

SMA(x, 20)

[1] NA NA NA NA NA NA NA NA NA NA

[11] NA NA NA NA NA NA NA NA NA 3.2420

[21] 3.2600 3.2750 3.2905 3.3070 3.3245 3.3360 3.3465 3.3585 3.3665 3.3765

[31] 3.3895 3.4035 3.4210 3.4330 3.4425 3.4500 3.4570 3.4610 3.4650 3.4685

[41] 3.4720 3.4790 3.4830 3.4860 3.4850 3.4855 3.4840 3.4800 3.4770 3.4705

mean(x[1:20])

[1] 3.242

mean(x[2:21])

[1] 3.26

Data Analysis using R Statistical Package

E-Content For the courses AEC - PGS 504 , Basic Laboratory Techniques , PAJANCOA & STAT 522, Data Analysis using R Statistical Packages , KAU

By

Mr.K.C.Ayyoob, Assistant Professor (Agricultural Statistics)KAU, Kerala-680656

Dr.S.Saravanan, Associate Professor (Computer Science)PAJANCOA & RI, Puducherry-609603

R programming language is a free and open-source software used for all kinds of data analysis, statistics, and visualization projects. These notes are an introduction R software for an introductory and advanced level statistics course.

R comes with several built-in data sets, which are generally used as demo data for demonstrating R functions. Majority of the examples here are shown using the built in data sets in R

R - Generating Frequency Tables

R - Summary

F Test

Correlation

ANOVA Without Interactions Between Factors (N, P & K)

ANOVA With Interactions Between Factors (N, P & K)

ANCOVA (Analyisis of Covariance)

Simple Linear Regression (lm - Linear Model)

Multiple Linear Regression (lm - Linear Model)

Maximum Likelihood Estimation (MLE)

Econometrics - Logit Model (Logistic Regression / Logit Regression), Probit Model

Stepwise Regression

Kolmogorov-Smirnov Test

Wilcoxon Signed Rank Test

Mann-Whitney U Test (Two Sample Unpaired Wilcoxon Test)

Kruskal-Wallis Test

McNemar's Test

Factor Analysis

Principal Component Analysis

R - Example For Creating / Testing Time Series Objects

R - Example For AutoRegressive (AR) Modelling of Time Series & Predicting

R - Example For MA (Moving Average)

Mr.K.C.Ayyoob, Assistant Professor (Agricultural Statistics)
KAU, Kerala-680656

Dr.S.Saravanan, Associate Professor (Computer Science)
PAJANCOA & RI, Puducherry-609603