Regression with a binary dependent variable

Author

Zahid Asghar

School of Economics, QAU, Islamabad

Regression with a binary dependent variable

So far the dependent variable (\(Y\)) has been continuous:

  • district-wide average test score
  • economic growth rate
  • earnings

What if \(Y\) is binary?

  • \(Y\) = get into college, or not; \(X\) = years of education

- \(Y\) = person smokes, or not; \(X\) = income

- \(Y\) = mortgage application is accepted, or not; \(X\) = income, house characteristics, marital status, race

Examples of Binary Dependent Variables
Topic Dummy Dependent Variable Description
Labour Force Participation Inlabourforce

0 if out of LF

1 if in the LF
Choice of Occupation Managerial

0 if not managerial

1 if managerial
Firm Location Shoppingmall

0 if not in the shopping mall

1 if in the shopping mall
Union Membership Union

0 if not a union member

1 if a union member
Retirement Retired

0 if not retired

1 if retired
Use of Seat Belts Seatbeltused

0 if does not use seat belt

1 if uses seat belt

Main concepts

We review the following concepts:

  • the linear probability model
  • the Probit model
  • the Logit model
  • maximum likelihood estimation of nonlinear regression models

Scatter plot of the Boston HMDA data

Code 
library(haven)
hmda <- read_stata("hmda_sw.dta")
# create the variables of interest
hmda$deny = ifelse(hmda$s7==3, 1, 0)
hmda$pi_ratio = hmda$s46/100
hmda$black = ifelse(hmda$s13==3, 1, 0)
library(ggplot2)
figure11_1 = ggplot(data=hmda, aes(x=pi_ratio, y=deny))
figure11_1 = figure11_1 +
geom_point(alpha=0.2, size=1) +
geom_smooth(method="lm", se=F) +
ylim(-0.4,1.4) + scale_x_continuous(breaks = seq(0,3,0.5)) +
labs(x="P/I Ratio", y="Deny") +
annotate("text",x=2.5,y=0.9,label="Mortgage denied") +
annotate("text",x=2.5,y=-0.1,label="Mortgage approved") +
theme_bw()
figure11_1
`geom_smooth()` using formula 'y ~ x'
Warning: Removed 15 rows containing missing values (geom_smooth).

Trimming the data

Code 
# sample of the HMDA data with 150 observations
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
Code 
hmda_s1 = dplyr::sample_n(filter(hmda,s13==3 & s7==3 & s46/100<1), 35)
hmda_s2 = dplyr::sample_n(filter(hmda,s13==3 & s46/100<1), 15)
hmda_s3 = dplyr::sample_n(filter(hmda,s13==5 & s46/100<1), 50)
hmda_s = rbind(hmda_s1, hmda_s2, hmda_s3)
# create the variables of interest
hmda_s$deny = ifelse(hmda_s$s7==3, 1, 0)
hmda_s$pi_ratio = hmda_s$s46/100
hmda_s$black = ifelse(hmda_s$s13==3, 1, 0)
figure11_1mod = ggplot(data=hmda_s, aes(x=pi_ratio, y=deny))
figure11_1mod = figure11_1mod +
geom_point(size=2, shape=1) +
geom_smooth(method="lm", se=F) +
ylim(-0.4,1.4) + scale_x_continuous(breaks = seq(0,1,0.1)) +
labs(x="P/I Ratio", y="Deny") +
annotate("text",x=0.6,y=1.1,label="Mortgage denied") +
annotate("text",x=0.6,y=-0.1,label="Mortgage approved") +
annotate("label",x=0.3,y=0.5,label="LPM") +
theme_bw()
figure11_1mod
`geom_smooth()` using formula 'y ~ x'

LPM

Equation of the LPM regression line: \(\begin{align} deny = \beta_0 + \beta_1 \times P/I\ ratio + u. \tag{11.1} \end{align}\)

\(\begin{align} \widehat{deny} = -\underset{(0.032)}{0.080} + \underset{(0.098)}{0.604} P/I \ ratio. \tag{11.2} \end{align}\)

Code 
library(lmtest)
Loading required package: zoo

Attaching package: 'zoo'
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
Code 
library(sandwich)
LPM = lm(deny ~ pi_ratio, data=hmda)
coeftest(LPM, vcov = sandwich)

t test of coefficients:

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.079910   0.031953 -2.5008  0.01246 *  
pi_ratio     0.603535   0.098441  6.1309 1.02e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LPM inclduing race variable

\(\begin{align} \widehat{deny} =& \, -\underset{(0.029)}{0.091} + \underset{(0.089)}{0.559} P/I \ ratio + \underset{(0.025)}{0.177} black. \tag{11.3} \end{align}\)

Code 
LPM2 = lm(deny ~ pi_ratio + black, data=hmda)
coeftest(LPM2, vcov=sandwich)

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -0.090514   0.028582 -3.1669  0.001561 ** 
pi_ratio     0.559195   0.088610  6.3107 3.303e-10 ***
black        0.177428   0.024931  7.1169 1.455e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Logit and Probit model

Code 
library(ggplot2)
ggplot(data.frame(x = c(-5,5)), aes(x=x)) + 
  stat_function(fun = pnorm, aes(colour = "Probit")) + 
  stat_function(fun = plogis, aes(colour = "Logit")) + 
  theme_bw() + 
  scale_colour_manual(name = "Function G",values = c("red", "blue")) +
  scale_y_continuous(name = "Pr(y = 1 | x)")

Probit Model, Predicted Probabilities and Estimated Effects

Assume that \(Y\) is a binary variable. The model \[Y= \beta_0 + \beta_1 + X_1 + \beta_2 X_2 + \dots + \beta_k X_k + u\] with \(P(Y = 1 \vert X_1, X_2, \dots ,X_k)\\ = \Phi(\beta_0 + \beta_1 + X_1 + \beta_2 X_2 + \dots + \beta_k X_k)\) is the population Probit model with multiple regressors \(X_1, X_2, \dots, X_k\) and \(\Phi(\cdot)\) is the cumulative standard normal distribution function.

Logit regression

The population Logit regression function is \[\begin{align*} P(Y=1\vert X_1, X_2, \dots, X_k) \\=& \, F(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k) \\ =& \, \frac{1}{1+e^{-(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k)}}. \end{align*}\] The idea is similar to Probit regression except that a different CDF is used: \(F(x) = \frac{1}{1+e^{-x}}\)

is the CDF of a standard logistically distributed random variable.

Probit regression

\(\begin{align} E(Y\vert X) = P(Y=1\vert X) = \Phi(\beta_0 + \beta_1 X). \tag{11.4} \end{align}\) \(\beta_0 + \beta_1 X\) in (11.4) plays a role of a quantile \(z\).Remember that \(\Phi(z) = P(Z \leq z) \ , \ Z \sim \mathcal{N}(0,1)\) such that the Probit coefficient \(\beta_1\) in (11.4) is the change in \(z\) associated with a one unit change in \(X\). Although the effect on \(z\) of a change in \(X\) is linear, the link between \(z\) and the dependent variable \(Y\) is nonlinear since \(\Phi\) is a nonlinear function of \(X\). .

Probit chart code

Code 
figure11_3mod = ggplot(data=hmda_s, aes(x=pi_ratio, y=deny))
figure11_3probit = figure11_3mod +
geom_point(shape=1, size=2) +
geom_smooth(method="glm",

method.args=list(family=binomial(link="probit")),
se=F) +
ylim(-0.4,1.4) +
scale_x_continuous(breaks = seq(0,1,0.1)) +
labs(x="P/I Ratio", y="Deny") +
annotate("label",x=0.3,y=0.5,label="Probit Regression") +
theme_bw()

Probit regression

\(\begin{align} \widehat{P(deny\vert P/I \ ratio, black)} = \Phi (-\underset{(0.18)}{2.26} + \underset{(0.50)}{2.74} P/I \ ratio + \underset{(0.08)}{0.71} black) \end{align}\)

Code 
probit = glm(
deny ~ black + pi_ratio,
data=hmda,
family=binomial(link="logit"))
coeftest(probit, vcov=sandwich)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
(Intercept) -4.12556    0.34576 -11.9320 < 2.2e-16 ***
black        1.27278    0.14607   8.7136 < 2.2e-16 ***
pi_ratio     5.37036    0.96315   5.5758 2.464e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

LPM , Logit, Probit

Code 
# rename and create some variables for regression
hmda$hse_inc = hmda$s45/100
hmda$loan_val = hmda$s6/hmda$s50
hmda$ccred = hmda$s43
hmda$mcred = hmda$s42
hmda$pubrec = ifelse(hmda$s44>0,1,0)
hmda$denpmi = ifelse(hmda$s53==1,1,0)
hmda$selfemp = ifelse(hmda$s27a==1,1,0)
hmda$married = ifelse(hmda$s23a=="M",1,0)
hmda$single = ifelse(hmda$married==0,1,0)
hmda$hischl = ifelse(hmda$school>=12,1,0)
hmda$probunmp = hmda$uria
hmda$condo = ifelse(hmda$s51 == 1,1,0)
Code 
# LPM
model_1 = lm(
deny ~ black + pi_ratio + hse_inc + loan_val + ccred +
mcred + pubrec + + denpmi + selfemp,
data=hmda)
# Logit
model_2 = glm(
deny ~ black + pi_ratio + hse_inc + loan_val + ccred +
mcred + pubrec + + denpmi + selfemp,
data=hmda,
family=binomial(link="logit"))
# Probit
model_3 = glm(
deny ~ black + pi_ratio + hse_inc + loan_val + ccred +
mcred + pubrec + + denpmi + selfemp,
data=hmda,
family=binomial(link="probit"))
Code 
models<-list(model_1,model_2,model_3)
library(huxtable)

Attaching package: 'huxtable'
The following object is masked from 'package:dplyr':

    add_rownames
The following object is masked from 'package:ggplot2':

    theme_grey
Code 
library(modelsummary)
modelsummary(models,estimate = "{estimate}{stars}", output="huxtable")
Model 1Model 2Model 3
(Intercept)-0.241***-6.830***-3.519***
(0.032)(0.541)(0.266)
black0.090***0.726***0.408***
(0.017)(0.175)(0.096)
pi_ratio0.461***4.743***2.472***
(0.087)(1.046)(0.549)
hse_inc-0.054-0.191-0.238
(0.096)(1.234)(0.653)
loan_val0.092**1.784***0.764**
(0.034)(0.498)(0.251)
ccred0.031***0.290***0.154***
(0.004)(0.040)(0.021)
mcred0.023*0.307*0.161*
(0.011)(0.139)(0.073)
pubrec0.201***1.221***0.701***
(0.023)(0.203)(0.117)
denpmi0.712***4.518***2.542***
(0.041)(0.550)(0.281)
selfemp0.060***0.672**0.354**
(0.018)(0.210)(0.111)
Num.Obs.238023802380
R20.258
R2 Adj.0.255
AIC711.71302.11305.7
BIC775.21359.81363.4
Log.Lik.-344.855-641.046-642.847
F91.47932.29835.294
RMSE0.280.280.28

mroz data from AER/Woolridge

Code 
data(mroz, package = "wooldridge")
plot(factor(inlf) ~ age, data = mroz, 
     ylevels = 2:1,
     ylab = "in labor force?")

Code 
LPM = lm(inlf ~ nwifeinc + educ + exper 
         + I(exper^2) + age +I(age^2) + kidslt6, mroz)
summary(LPM)

Call:
lm(formula = inlf ~ nwifeinc + educ + exper + I(exper^2) + age + 
    I(age^2) + kidslt6, data = mroz)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.94194 -0.37773  0.08935  0.34283  0.97979 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.3219554  0.4863667   0.662  0.50820    
nwifeinc    -0.0034271  0.0014531  -2.358  0.01861 *  
educ         0.0374662  0.0073476   5.099 4.33e-07 ***
exper        0.0382568  0.0057700   6.630 6.44e-11 ***
I(exper^2)  -0.0005649  0.0001895  -2.981  0.00296 ** 
age         -0.0011177  0.0225100  -0.050  0.96041    
I(age^2)    -0.0001822  0.0002581  -0.706  0.48044    
kidslt6     -0.2603675  0.0340826  -7.639 6.72e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4273 on 745 degrees of freedom
Multiple R-squared:  0.2637,    Adjusted R-squared:  0.2568 
F-statistic: 38.13 on 7 and 745 DF,  p-value: < 2.2e-16

Code 
pr = predict(LPM)
plot(pr[order(pr)],ylab = "p(inlf = 1)")
abline(a = 0, b = 0, col = "red")
abline(a = 1, b = 0, col = "red")

Code 
library(dplyr)
library(ggplot2)
mroz %<>% 
  # classify age into 3 and huswage into 2 classes
  mutate(age_fct = cut(age,breaks = 3,labels = FALSE),
         huswage_fct = cut(huswage, breaks = 2,labels = FALSE)) %>%
  mutate(classes = paste0("age_",age_fct,"_hus_",huswage_fct))

LPM_saturated = mroz %>%
  lm(inlf ~ classes, data = .)

mroz$pred <- predict(LPM_saturated)

ggplot(mroz[order(mroz$pred),], aes(x = 1:nrow(mroz),y = pred,color = classes)) + 
  geom_point() + 
  theme_bw() + 
  scale_y_continuous(limits = c(0,1), name = "p(inlf)") +
  ggtitle("LPM in a Saturated Model is Perfectly Fine")

Interpretation of coefficients

Code 
probit <- glm(inlf ~ age, 
                    data = mroz, 
                    family = binomial(link = "probit"))

logit <- glm(inlf ~ age, 
                    data = mroz, 
                    family = binomial(link = "logit"))
modelsummary::modelsummary(list("probit" = probit,"logit" = logit))
probit logit
(Intercept) 0.707 1.136
(0.248) (0.398)
age −0.013 −0.020
(0.006) (0.009)
Num.Obs. 753 753
AIC 1028.9 1028.9
BIC 1038.1 1038.1
Log.Lik. −512.442 −512.431
F 4.828 4.858
RMSE 0.49 0.49

shinyapp

Code 
library(mfx)
Loading required package: MASS

Attaching package: 'MASS'
The following object is masked from 'package:dplyr':

    select
Loading required package: betareg
Code 
f <- "inlf ~ age + kidslt6 + nwifeinc" # setup a formula
glms <- list()
glms$probit <- glm(formula = f, 
                    data = mroz, 
                    family = binomial(link = "probit"))
glms$logit <- glm(formula = f, 
                    data = mroz, 
                    family = binomial(link = "logit"))
# now the marginal effects versions
glms$probitMean <- mfx::probitmfx(formula = f, 
                    data = mroz, atmean = TRUE)
glms$probitAvg <- mfx::probitmfx(formula = f, 
                    data = mroz, atmean = FALSE)
glms$logitMean <- mfx::logitmfx(formula = f, 
                    data = mroz, atmean = TRUE)
glms$logitAvg <- mfx::logitmfx(formula = f, 
                    data = mroz, atmean = FALSE)
modelsummary::modelsummary(glms, 
                           stars = TRUE,
                           gof_omit = "AIC|BIC",
                           title = "Logit and Probit estimates and marginal effects evaluated at mean of x or as sample average of effects")
Warning: 
There are duplicate term names in the table.

The `shape` argument of the `modelsummary` function can be used to print
related terms together. The `group_map` argument can be used to reorder,
subset, and rename group identifiers. See `?modelsummary` for details.

You can find the group identifier to use in the `shape` argument by calling
`get_estimates()` on one of your models. Candidates include: group, component
Logit and Probit estimates and marginal effects evaluated at mean of x or as sample average of effects
probit logit probitMean probitAvg logitMean logitAvg
(Intercept) 2.080*** 3.394*** 2.080*** 2.080*** 3.394*** 3.394***
(0.309) (0.516) (0.309) (0.309) (0.516) (0.516)
age −0.035*** −0.057*** −0.014*** −0.013*** −0.014*** −0.013***
−0.035*** −0.057*** −0.014*** −0.013*** −0.014*** −0.057***
−0.035*** −0.057*** −0.014*** −0.013*** −0.057*** −0.013***
−0.035*** −0.057*** −0.014*** −0.013*** −0.057*** −0.057***
−0.035*** −0.057*** −0.014*** −0.035*** −0.014*** −0.013***
−0.035*** −0.057*** −0.014*** −0.035*** −0.014*** −0.057***
−0.035*** −0.057*** −0.014*** −0.035*** −0.057*** −0.013***
−0.035*** −0.057*** −0.014*** −0.035*** −0.057*** −0.057***
−0.035*** −0.057*** −0.035*** −0.013*** −0.014*** −0.013***
−0.035*** −0.057*** −0.035*** −0.013*** −0.014*** −0.057***
−0.035*** −0.057*** −0.035*** −0.013*** −0.057*** −0.013***
−0.035*** −0.057*** −0.035*** −0.013*** −0.057*** −0.057***
−0.035*** −0.057*** −0.035*** −0.035*** −0.014*** −0.013***
−0.035*** −0.057*** −0.035*** −0.035*** −0.014*** −0.057***
−0.035*** −0.057*** −0.035*** −0.035*** −0.057*** −0.013***
−0.035*** −0.057*** −0.035*** −0.035*** −0.057*** −0.057***
(0.007) (0.011) (0.003) (0.002) (0.003) (0.003)
(0.007) (0.011) (0.003) (0.002) (0.003) (0.011)
(0.007) (0.011) (0.003) (0.002) (0.011) (0.003)
(0.007) (0.011) (0.003) (0.002) (0.011) (0.011)
(0.007) (0.011) (0.003) (0.007) (0.003) (0.003)
(0.007) (0.011) (0.003) (0.007) (0.003) (0.011)
(0.007) (0.011) (0.003) (0.007) (0.011) (0.003)
(0.007) (0.011) (0.003) (0.007) (0.011) (0.011)
(0.007) (0.011) (0.007) (0.002) (0.003) (0.003)
(0.007) (0.011) (0.007) (0.002) (0.003) (0.011)
(0.007) (0.011) (0.007) (0.002) (0.011) (0.003)
(0.007) (0.011) (0.007) (0.002) (0.011) (0.011)
(0.007) (0.011) (0.007) (0.007) (0.003) (0.003)
(0.007) (0.011) (0.007) (0.007) (0.003) (0.011)
(0.007) (0.011) (0.007) (0.007) (0.011) (0.003)
(0.007) (0.011) (0.007) (0.007) (0.011) (0.011)
kidslt6 −0.800*** −1.313*** −0.314*** −0.290*** −0.322*** −0.292***
−0.800*** −1.313*** −0.314*** −0.290*** −0.322*** −1.313***
−0.800*** −1.313*** −0.314*** −0.290*** −1.313*** −0.292***
−0.800*** −1.313*** −0.314*** −0.290*** −1.313*** −1.313***
−0.800*** −1.313*** −0.314*** −0.800*** −0.322*** −0.292***
−0.800*** −1.313*** −0.314*** −0.800*** −0.322*** −1.313***
−0.800*** −1.313*** −0.314*** −0.800*** −1.313*** −0.292***
−0.800*** −1.313*** −0.314*** −0.800*** −1.313*** −1.313***
−0.800*** −1.313*** −0.800*** −0.290*** −0.322*** −0.292***
−0.800*** −1.313*** −0.800*** −0.290*** −0.322*** −1.313***
−0.800*** −1.313*** −0.800*** −0.290*** −1.313*** −0.292***
−0.800*** −1.313*** −0.800*** −0.290*** −1.313*** −1.313***
−0.800*** −1.313*** −0.800*** −0.800*** −0.322*** −0.292***
−0.800*** −1.313*** −0.800*** −0.800*** −0.322*** −1.313***
−0.800*** −1.313*** −0.800*** −0.800*** −1.313*** −0.292***
−0.800*** −1.313*** −0.800*** −0.800*** −1.313*** −1.313***
(0.111) (0.188) (0.044) (0.036) (0.046) (0.047)
(0.111) (0.188) (0.044) (0.036) (0.046) (0.188)
(0.111) (0.188) (0.044) (0.036) (0.188) (0.047)
(0.111) (0.188) (0.044) (0.036) (0.188) (0.188)
(0.111) (0.188) (0.044) (0.111) (0.046) (0.047)
(0.111) (0.188) (0.044) (0.111) (0.046) (0.188)
(0.111) (0.188) (0.044) (0.111) (0.188) (0.047)
(0.111) (0.188) (0.044) (0.111) (0.188) (0.188)
(0.111) (0.188) (0.111) (0.036) (0.046) (0.047)
(0.111) (0.188) (0.111) (0.036) (0.046) (0.188)
(0.111) (0.188) (0.111) (0.036) (0.188) (0.047)
(0.111) (0.188) (0.111) (0.036) (0.188) (0.188)
(0.111) (0.188) (0.111) (0.111) (0.046) (0.047)
(0.111) (0.188) (0.111) (0.111) (0.046) (0.188)
(0.111) (0.188) (0.111) (0.111) (0.188) (0.047)
(0.111) (0.188) (0.111) (0.111) (0.188) (0.188)
nwifeinc −0.011** −0.019** −0.004** −0.004** −0.005** −0.004**
−0.011** −0.019** −0.004** −0.004** −0.005** −0.019**
−0.011** −0.019** −0.004** −0.004** −0.019** −0.004**
−0.011** −0.019** −0.004** −0.004** −0.019** −0.019**
−0.011** −0.019** −0.004** −0.011** −0.005** −0.004**
−0.011** −0.019** −0.004** −0.011** −0.005** −0.019**
−0.011** −0.019** −0.004** −0.011** −0.019** −0.004**
−0.011** −0.019** −0.004** −0.011** −0.019** −0.019**
−0.011** −0.019** −0.011** −0.004** −0.005** −0.004**
−0.011** −0.019** −0.011** −0.004** −0.005** −0.019**
−0.011** −0.019** −0.011** −0.004** −0.019** −0.004**
−0.011** −0.019** −0.011** −0.004** −0.019** −0.019**
−0.011** −0.019** −0.011** −0.011** −0.005** −0.004**
−0.011** −0.019** −0.011** −0.011** −0.005** −0.019**
−0.011** −0.019** −0.011** −0.011** −0.019** −0.004**
−0.011** −0.019** −0.011** −0.011** −0.019** −0.019**
(0.004) (0.007) (0.002) (0.001) (0.002) (0.002)
(0.004) (0.007) (0.002) (0.001) (0.002) (0.007)
(0.004) (0.007) (0.002) (0.001) (0.007) (0.002)
(0.004) (0.007) (0.002) (0.001) (0.007) (0.007)
(0.004) (0.007) (0.002) (0.004) (0.002) (0.002)
(0.004) (0.007) (0.002) (0.004) (0.002) (0.007)
(0.004) (0.007) (0.002) (0.004) (0.007) (0.002)
(0.004) (0.007) (0.002) (0.004) (0.007) (0.007)
(0.004) (0.007) (0.004) (0.001) (0.002) (0.002)
(0.004) (0.007) (0.004) (0.001) (0.002) (0.007)
(0.004) (0.007) (0.004) (0.001) (0.007) (0.002)
(0.004) (0.007) (0.004) (0.001) (0.007) (0.007)
(0.004) (0.007) (0.004) (0.004) (0.002) (0.002)
(0.004) (0.007) (0.004) (0.004) (0.002) (0.007)
(0.004) (0.007) (0.004) (0.004) (0.007) (0.002)
(0.004) (0.007) (0.004) (0.004) (0.007) (0.007)
Num.Obs. 753 753 753 753 753 753
Log.Lik. −478.395 −478.377
F 21.784 20.280
RMSE 0.47 0.47 0.47 0.47 0.47 0.47
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

This is why R

Code 
library(leaflet)
leaflet() %>%
  addTiles() %>%
  addMarkers(lng = 73.136946, lat =33.748294 ,
             popup = "School of Economics, QAU, Islamabad")