lect5

linear regression

table of content of lecture-5

• author’s details
  • author’s image
  • author’s website
• setting R environment
• multivariate linear regression analysis
  • generating dummy data
  • performing regression analysis
  • generating dummy data for prediction
  • predicting output using regression model
• employee data analysis
  • loading required function
  • load employee data
  • data visualisation
  • removing data with NA
  • calculating age & recoding gender
  • linear regression for employee’s data
  • data generation for prediction
  • predicting data using regression model
• closing R environment

author’s details

author’s image

author’s website

https://kiritved.com

setting R environment

rm(list=ls(all.names = T))
setwd("d:/met/met_lect5")
if(! require("pacman")) install.packages("pacman")

Loading required package: pacman

library("pacman")
p_load("tidyverse")
p_loaded()

 [1] "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
 [7] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "pacman"   

multivariate linear regression analysis

generating dummy data

n=10
  x1=runif(n,-5,5)
  x2=rnorm(n,50,10)
  x3=rnorm(n,20,5)
y=23.67*x1+45.87*x2+92.56*x3+runif(n,-50,50)
dd=data.frame(x1,x2,x3,y)
  
  
dd$y=y
dd

x1	x2	x3	y
-2.1104386	49.02694	19.65088	4042.283
-2.9485693	55.26326	11.91119	3590.408
-4.7167570	40.75296	20.60386	3661.059
-3.7579141	47.60471	21.01700	4080.431
2.8434707	35.68743	13.10243	2907.542
-2.3514930	50.16462	24.34514	4465.004
-2.5106343	46.30577	18.88831	3835.799
0.4568921	35.02861	22.03809	3687.804
-2.7049528	39.31916	18.37292	3445.016
1.4721070	62.88295	18.95012	4643.266

performing regression analysis

mylm=lm(y~.,dd)
mylm;summary(mylm)

Call:
lm(formula = y ~ ., data = dd)

Coefficients:
(Intercept)           x1           x2           x3  
      73.94        19.06        44.95        90.86  

Call:
lm(formula = y ~ ., data = dd)

Residuals:
     Min       1Q   Median       3Q      Max 
-31.0480 -14.9752  -0.4211  17.4881  28.6511 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   73.945     68.044   1.087  0.31888    
x1            19.060      3.892   4.898  0.00272 ** 
x2            44.950      1.024  43.899 9.36e-09 ***
x3            90.861      2.468  36.821 2.68e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 26.92 on 6 degrees of freedom
Multiple R-squared:  0.9981,    Adjusted R-squared:  0.9971 
F-statistic:  1046 on 3 and 6 DF,  p-value: 1.52e-08

generating dummy data for prediction

n=3
x1=runif(n,min(dd$x1),max(dd$x1))
x2=runif(n,min(dd$x2),max(dd$x2))
x3=runif(n,min(dd$x3),max(dd$x3))
ddp=data.frame(x1,x2,x3)
ddp

x1	x2	x3
2.267781	47.78700	12.19451
-3.972028	61.50773	15.33674
-3.285195	51.63996	17.15774

predicting output using regression model

mypred=predict(mylm,ddp)
#mypred
ddp$y=mypred
ddp

x1	x2	x3	y
2.267781	47.78700	12.19451	3373.201
-3.972028	61.50773	15.33674	4156.526
-3.285195	51.63996	17.15774	3891.518

employee data analysis

loading required function

# original order
#id,gender,bdate,educ,jobcat,salbegin,jobtime,prevexp,minority,salary
f1=function(dt,n=30){
  ndt=as.Date(dt, "%Y") + years(n)
  return (ndt)
}
f2=function(dt){
  rd=Sys.Date()-dt
  rd=rd/365.25
  rd=as.numeric(rd)
  rd=ceiling(rd)
  return (rd)
}

load employee data

d=read.csv("employee data.csv")
d$bdate=mdy(d$bdate)

Warning: 1 failed to parse.

head(d)

id	gender	bdate	educ	jobcat	salbegin	jobtime	prevexp	minority	salary
1	m	1952-02-03	15	3	27000	98	144	0	57000
2	m	1958-05-23	16	1	18750	98	36	0	40200
3	f	1929-07-26	12	1	12000	98	381	0	21450
4	f	1947-04-15	8	1	13200	98	190	0	21900
5	m	1955-02-09	15	1	21000	98	138	0	45000
6	m	1958-08-22	15	1	13500	98	67	0	32100

tail(d)

	id	gender	bdate	educ	jobcat	salbegin	jobtime	prevexp	minority	salary
469	469	f	1964-06-01	15	1	13950	64	57	0	25200
470	470	m	1964-01-22	12	1	15750	64	69	1	26250
471	471	m	1966-08-03	15	1	15750	64	32	1	26400
472	472	m	1966-02-21	15	1	15750	63	46	0	39150
473	473	f	1937-11-25	12	1	12750	63	139	0	21450
474	474	f	1968-11-05	12	1	14250	63	9	0	29400

str(d)

'data.frame':   474 obs. of  10 variables:
 $ id      : int  1 2 3 4 5 6 7 8 9 10 ...
 $ gender  : chr  "m" "m" "f" "f" ...
 $ bdate   : Date, format: "1952-02-03" "1958-05-23" ...
 $ educ    : int  15 16 12 8 15 15 15 12 15 12 ...
 $ jobcat  : int  3 1 1 1 1 1 1 1 1 1 ...
 $ salbegin: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
 $ jobtime : int  98 98 98 98 98 98 98 98 98 98 ...
 $ prevexp : int  144 36 381 190 138 67 114 0 115 244 ...
 $ minority: int  0 0 0 0 0 0 0 0 0 0 ...
 $ salary  : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...

summary(d)

       id           gender              bdate                 educ      
 Min.   :  1.0   Length:474         Min.   :1929-02-10   Min.   : 8.00  
 1st Qu.:119.2   Class :character   1st Qu.:1948-01-03   1st Qu.:12.00  
 Median :237.5   Mode  :character   Median :1962-01-23   Median :12.00  
 Mean   :237.5                      Mean   :1956-10-08   Mean   :13.49  
 3rd Qu.:355.8                      3rd Qu.:1965-07-06   3rd Qu.:15.00  
 Max.   :474.0                      Max.   :1971-02-10   Max.   :21.00  
                                    NA's   :1                           
     jobcat         salbegin        jobtime         prevexp      
 Min.   :1.000   Min.   : 9000   Min.   :63.00   Min.   :  0.00  
 1st Qu.:1.000   1st Qu.:12488   1st Qu.:72.00   1st Qu.: 19.25  
 Median :1.000   Median :15000   Median :81.00   Median : 55.00  
 Mean   :1.411   Mean   :17016   Mean   :81.11   Mean   : 95.86  
 3rd Qu.:1.000   3rd Qu.:17490   3rd Qu.:90.00   3rd Qu.:138.75  
 Max.   :3.000   Max.   :79980   Max.   :98.00   Max.   :476.00  
                                                                 
    minority          salary      
 Min.   :0.0000   Min.   : 15750  
 1st Qu.:0.0000   1st Qu.: 24000  
 Median :0.0000   Median : 28875  
 Mean   :0.2194   Mean   : 34420  
 3rd Qu.:0.0000   3rd Qu.: 36938  
 Max.   :1.0000   Max.   :135000  
                                  

data visualisation

## number of male & female employee along with %
mybarplot=function(df,cn,ttl="bar plot",xlbl,ylbl){
  colnm=colnames(df)[cn]
  v=df|>select(c(cn))
  t=table(v)
  ttl=paste(ttl,"for",colnm,sep=" ")
  b=barplot(t,col="lightgreen",main=ttl,xlab=colnm)
  text(b,min(t)*1.1,labels=as.character(t),col="red")
  #return (T)
}
mypie=function(df,cn,ttl="pie chart"){
  colnm=colnames(df)[cn]
  v=df|>select(c(cn))
  
  t=table(v)

  revnm=paste(names(t),t,sep="--")
  ttl=paste(ttl,"for",colnm,sep=" ")
    pie(t,labels=revnm,main=ttl)
}
dtmp=d|>group_by(gender)|>summarise(n=n())|>mutate(nper=100*n/sum(n))
dtmp

gender	n	nper
f	216	45.56962
m	258	54.43038

mybarplot(d,2)

Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
  # Was:
  data %>% select(cn)

  # Now:
  data %>% select(all_of(cn))

See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.

mybarplot(d,4)

mybarplot(d,5)

mybarplot(d,9)

mypie(d,2)

mypie(d,4)

mypie(d,5)

mypie(d,9)

plot(d$salbegin,d$salary,col="lightgreen")

boxplot(salary~gender,d,col="lightgreen",main="boxplot of salary vs gender")

boxplot(salary~jobcat,d,col="lightgreen",main="boxplot of salary vs jobcategory")

boxplot(salary~educ,d,col="lightgreen",main="boxplot of salary vs educ")

boxplot(salary~minority,d,col="lightgreen",main="boxplot of salary vs minority")

removing data with NA

tmp=complete.cases(d)
d=d[tmp==T,]
summary(d)

       id           gender              bdate                 educ      
 Min.   :  1.0   Length:473         Min.   :1929-02-10   Min.   : 8.00  
 1st Qu.:119.0   Class :character   1st Qu.:1948-01-03   1st Qu.:12.00  
 Median :237.0   Mode  :character   Median :1962-01-23   Median :12.00  
 Mean   :237.1                      Mean   :1956-10-08   Mean   :13.49  
 3rd Qu.:355.0                      3rd Qu.:1965-07-06   3rd Qu.:15.00  
 Max.   :474.0                      Max.   :1971-02-10   Max.   :21.00  
     jobcat         salbegin        jobtime         prevexp      
 Min.   :1.000   Min.   : 9000   Min.   :63.00   Min.   :  0.00  
 1st Qu.:1.000   1st Qu.:12450   1st Qu.:72.00   1st Qu.: 19.00  
 Median :1.000   Median :15000   Median :81.00   Median : 55.00  
 Mean   :1.412   Mean   :17009   Mean   :81.14   Mean   : 95.95  
 3rd Qu.:1.000   3rd Qu.:17490   3rd Qu.:90.00   3rd Qu.:139.00  
 Max.   :3.000   Max.   :79980   Max.   :98.00   Max.   :476.00  
    minority          salary      
 Min.   :0.0000   Min.   : 15750  
 1st Qu.:0.0000   1st Qu.: 24000  
 Median :0.0000   Median : 28800  
 Mean   :0.2199   Mean   : 34418  
 3rd Qu.:0.0000   3rd Qu.: 37050  
 Max.   :1.0000   Max.   :135000  

calculating age & recoding gender

d1=d|> mutate(
  m1=f1(bdate,30),
  age=f2(m1),
  gn=recode(gender,"m"=1,"f"=0)
  )|> select(-c(1,2,3,11))|>select(c(1:6,8:9),7)
#d1$bdate=f1(d1$bdate,30)
#write.csv(d1,"mycleaned.csv")
#d1$age=f2(d1$bdate)
#d1$gn=recode(d1$gender,"m"=1,"f"=0)
head(d1)

educ	jobcat	salbegin	jobtime	prevexp	minority	age	gn	salary
15	3	27000	98	144	0	42	1	57000
16	1	18750	98	36	0	36	1	40200
12	1	12000	98	381	0	65	0	21450
8	1	13200	98	190	0	47	0	21900
15	1	21000	98	138	0	39	1	45000
15	1	13500	98	67	0	36	1	32100

str(d1)

'data.frame':   473 obs. of  9 variables:
 $ educ    : int  15 16 12 8 15 15 15 12 15 12 ...
 $ jobcat  : int  3 1 1 1 1 1 1 1 1 1 ...
 $ salbegin: int  27000 18750 12000 13200 21000 13500 18750 9750 12750 13500 ...
 $ jobtime : int  98 98 98 98 98 98 98 98 98 98 ...
 $ prevexp : int  144 36 381 190 138 67 114 0 115 244 ...
 $ minority: int  0 0 0 0 0 0 0 0 0 0 ...
 $ age     : num  42 36 65 47 39 36 38 28 48 48 ...
 $ gn      : num  1 1 0 0 1 1 1 0 0 0 ...
 $ salary  : int  57000 40200 21450 21900 45000 32100 36000 21900 27900 24000 ...

linear regression for employee’s data

#head(d1)
#str(d1)
mylm=lm(salary~.,d1)
mylm

Call:
lm(formula = salary ~ ., data = d1)

Coefficients:
(Intercept)         educ       jobcat     salbegin      jobtime      prevexp  
 -11780.031      455.332     5753.745        1.329      154.809      -14.604  
   minority          age           gn  
   -970.473      -68.997     1804.891  

summary(mylm)

Call:
lm(formula = salary ~ ., data = d1)

Residuals:
   Min     1Q Median     3Q    Max 
-23214  -3018   -733   2655  46335 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.178e+04  3.252e+03  -3.622 0.000324 ***
educ         4.553e+02  1.541e+02   2.956 0.003279 ** 
jobcat       5.754e+03  6.226e+02   9.242  < 2e-16 ***
salbegin     1.329e+00  7.037e-02  18.883  < 2e-16 ***
jobtime      1.548e+02  3.164e+01   4.893 1.37e-06 ***
prevexp     -1.460e+01  5.507e+00  -2.652 0.008278 ** 
minority    -9.705e+02  7.844e+02  -1.237 0.216659    
age         -6.900e+01  4.772e+01  -1.446 0.148892    
gn           1.805e+03  7.744e+02   2.331 0.020194 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6803 on 464 degrees of freedom
Multiple R-squared:  0.8443,    Adjusted R-squared:  0.8416 
F-statistic: 314.5 on 8 and 464 DF,  p-value: < 2.2e-16

data generation for prediction

cn=colnames(d)
n=10
id=1:n
gender=sample(c("m","f"),n,replace = T)
bdate=runif(n,min(d$bdate),max(d$bdate))
educ=sample(min(d$educ):max(d$educ),n,replace=T)
jobcat=sample(c(1:3),n,replace = T)
salbegin=runif(n,min(d$salbegin),max(d$salbegin))
salbegin=ceiling(salbegin)
prevexp=runif(n,min(d$jobtime),max(d$jobtime))
prevexp=ceiling(prevexp)
jobtime=runif(n,min(d$jobtime),max(d$jobtime))
jobtime=ceiling(jobtime)
minority=sample(c(0,1),n,replace = T)
dp=data.frame(id,gender,bdate=as.Date(bdate),educ,jobcat,salbegin,jobtime,prevexp,minority)
dp

id	gender	bdate	educ	jobcat	salbegin	jobtime	prevexp	minority
1	m	1945-12-16	11	1	27442	96	88	0
2	f	1944-02-03	21	3	46762	72	88	1
3	f	1970-04-16	20	2	31907	84	90	1
4	f	1941-10-07	12	3	19337	68	66	0
5	m	1962-12-15	21	2	74095	74	67	0
6	f	1956-01-24	11	2	22504	67	92	1
7	f	1954-06-16	9	3	78082	86	70	1
8	f	1969-10-02	8	1	23489	90	67	1
9	m	1943-07-16	17	2	14939	73	93	1
10	f	1933-08-28	17	3	11143	97	93	1

dp1=dp|> mutate(
 # m1=f1(bdate,30),
  age=f2(bdate),
  gn=recode(gender,"m"=1,"f"=0))  |> select(-c(1,2,3))
#d1$bdate=f1(d1$bdate,30)
#write.csv(d1,"mycleaned.csv")
#d1$age=f2(d1$bdate)
#d1$gn=recode(d1$gender,"m"=1,"f"=0)
head(dp1)

educ	jobcat	salbegin	jobtime	prevexp	minority	age	gn
11	1	27442	96	88	0	78	1
21	3	46762	72	88	1	80	0
20	2	31907	84	90	1	54	0
12	3	19337	68	66	0	83	0
21	2	74095	74	67	0	61	1
11	2	22504	67	92	1	68	0

str(dp1)

'data.frame':   10 obs. of  8 variables:
 $ educ    : int  11 21 20 12 21 11 9 8 17 17
 $ jobcat  : int  1 3 2 3 2 2 3 1 2 3
 $ salbegin: num  27442 46762 31907 19337 74095 ...
 $ jobtime : num  96 72 84 68 74 67 86 90 73 97
 $ prevexp : num  88 88 90 66 67 92 70 67 93 93
 $ minority: num  0 1 1 0 0 1 1 1 1 1
 $ age     : num  78 80 54 83 61 68 70 55 81 91
 $ gn      : num  1 0 0 0 1 0 0 0 1 0

predicting data using regression model

myp=predict(mylm,dp1)
#myp
dp1$salarypred=myp
dp1

educ	jobcat	salbegin	jobtime	prevexp	minority	age	gn	salarypred
11	1	27442	96	88	0	78	1	45446.53
21	3	46762	72	88	1	80	0	80550.64
20	2	31907	84	90	1	54	0	58224.92
12	3	19337	68	66	0	83	0	40476.31
21	2	74095	74	67	0	61	1	115819.11
11	2	22504	67	92	1	68	0	38005.46
9	3	78082	86	70	1	70	0	119824.29
8	1	23489	90	67	1	55	0	37017.24
17	2	14939	73	93	1	81	1	32507.40
17	3	11143	97	93	1	91	0	34437.65

closing R environment

p_unload(all)

The following packages have been unloaded:
lubridate, forcats, stringr, dplyr, purrr, readr, tidyr, tibble, ggplot2, tidyverse, pacman

rm(list=ls(all.names = T))