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

Code

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

Loading required package: pacman

Code

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

Code

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
0.1651063	56.11450	17.81261	4197.918
1.4214578	73.85352	20.48178	5357.737
4.8131084	38.79473	22.16092	3969.592
-2.7655074	52.01693	26.60316	4801.017
-4.3147171	39.34803	13.76159	3003.366
-1.7348955	45.78057	24.55680	4295.714
2.4674704	52.88396	14.64843	3883.709
3.9479441	40.81402	12.38804	3065.226
1.3113165	67.68993	33.45296	6214.195
2.1776584	55.64402	10.93243	3615.444

performing regression analysis

Code

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

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

Coefficients:
(Intercept)           x1           x2           x3  
     -15.39        22.48        46.57        91.65  

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

Residuals:
   Min     1Q Median     3Q    Max 
-44.19 -28.87   3.32  24.05  39.11 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -15.385     61.216  -0.251  0.80995    
x1            22.482      4.563   4.927  0.00264 ** 
x2            46.570      1.202  38.739 1.98e-08 ***
x3            91.647      2.005  45.720 7.33e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 39.14 on 6 degrees of freedom
Multiple R-squared:  0.999, Adjusted R-squared:  0.9985 
F-statistic:  1949 on 3 and 6 DF,  p-value: 2.357e-09

generating dummy data for prediction

Code

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
3.6928321	52.37799	18.71128
1.1813351	41.98070	15.95708
-0.5724092	51.26017	24.34673

predicting output using regression model

Code

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

x1	x2	x3	y
3.6928321	52.37799	18.71128	4221.719
1.1813351	41.98070	15.95708	3428.638
-0.5724092	51.26017	24.34673	4590.244

employee data analysis

loading required function

Code

# 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

Code

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

Warning: 1 failed to parse.

Code

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

Code

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

Code

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 ...

Code

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

Code

## 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

Code

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>.

Code

mybarplot(d,4)

Code

mybarplot(d,5)

Code

mybarplot(d,9)

Code

mypie(d,2)

Code

mypie(d,4)

Code

mypie(d,5)

Code

mypie(d,9)

Code

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

Code

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

Code

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

Code

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

Code

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

removing data with NA

Code

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

Code

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

Code

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

Code

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

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

Coefficients:
(Intercept)         educ       jobcat     salbegin      jobtime      prevexp  
 -11790.064      455.359     5753.970        1.329      154.774      -14.638  
   minority          age           gn  
   -970.220      -68.605     1806.708  

Code

summary(mylm)

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

Residuals:
   Min     1Q Median     3Q    Max 
-23213  -3026   -730   2658  46336 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.179e+04  3.252e+03  -3.625 0.000321 ***
educ         4.554e+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.882  < 2e-16 ***
jobtime      1.548e+02  3.164e+01   4.891 1.38e-06 ***
prevexp     -1.464e+01  5.507e+00  -2.658 0.008134 ** 
minority    -9.702e+02  7.845e+02  -1.237 0.216793    
age         -6.860e+01  4.772e+01  -1.438 0.151169    
gn           1.807e+03  7.744e+02   2.333 0.020077 *  
---
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.4 on 8 and 464 DF,  p-value: < 2.2e-16

data generation for prediction

Code

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	1957-03-18	11	2	23266	69	92	1
2	f	1951-02-22	11	1	50326	78	71	1
3	f	1953-03-03	9	2	78315	73	71	1
4	f	1962-09-04	21	1	32818	81	79	1
5	m	1958-04-29	19	2	61071	71	66	1
6	f	1946-11-01	19	2	9901	96	90	0
7	m	1935-12-29	16	2	17211	87	86	1
8	f	1963-09-22	9	1	22547	88	96	0
9	m	1963-10-13	10	2	76373	91	71	1
10	m	1970-05-07	19	2	41725	82	75	1

Code

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	2	23266	69	92	1	67	1
11	1	50326	78	71	1	73	0
9	2	78315	73	71	1	71	0
21	1	32818	81	79	1	62	0
19	2	61071	71	66	1	66	1
19	2	9901	96	90	0	78	0

Code

str(dp1)

'data.frame':   10 obs. of  8 variables:
 $ educ    : int  11 11 9 21 19 19 16 9 10 19
 $ jobcat  : int  2 1 2 1 2 2 2 1 2 2
 $ salbegin: num  23266 50326 78315 32818 61071 ...
 $ jobtime : num  69 78 73 81 71 96 87 88 91 82
 $ prevexp : num  92 71 71 79 66 90 86 96 71 75
 $ minority: num  1 1 1 1 1 0 1 0 1 1
 $ age     : num  67 73 71 62 66 78 88 61 61 54
 $ gn      : num  1 0 0 0 1 0 1 0 1 1

predicting data using regression model

Code

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

educ	jobcat	salbegin	jobtime	prevexp	minority	age	gn	salarypred
11	2	23266	69	92	1	67	1	41214.15
11	1	50326	78	71	1	73	0	70898.14
9	2	78315	73	71	1	71	0	112295.06
21	1	32818	81	79	1	62	0	53289.88
19	2	61071	71	66	1	66	1	95849.09
19	2	9901	96	90	0	78	0	29715.35
16	2	17211	87	86	1	88	1	36878.44
9	1	22547	88	96	0	61	0	36051.39
10	2	76373	91	71	1	61	1	115448.69
19	2	41725	82	75	1	54	1	72537.17

closing R environment

Code

p_unload(all)

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

Code

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