R Intermediate - Class 4
Tung Linh Ha
2024-05-27
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-library(Publish)## Loading required package: prodlimlibrary(rms)## Loading required package: Hmisc##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:table1':
##
## label, label<-, units## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(BMA)## Loading required package: survival## Loading required package: leaps## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.7-5)Task 1: Determine whether the chance to get a high salary differed between male and female professors
Task 1.1: Study Design:
- The cross-sectional investigation of 397 professors to determine whether the chance to get a high salary differed between male and female professors, accounting for potential confounding effects.
- Null hypothesis: Professors’ sexes were not independently associated with the chance to get a high salary.
- Alternative hypothesis: Professors’ sexes were independently associated with the chance to get a high salary
Task 1.2: Import data set
salary = read.csv("C:\\Users\\DELL\\OneDrive\\ANDREAS\\ACADEMICS_UNIVERSITY\\Year_4_2024\\Autumn\\32931\\Classes\\R_Int\\Professorial Salaries.csv")Task 1.3: Create a new variable ‘high.salary’ using the threshold of $120,000
salary$high.salary = ifelse(salary$Salary> 120000, 1, 0)
salary$Prof[salary$Rank == "Prof"] = "Professor"
salary$Prof[salary$Rank != "Prof"] = "AssProfessor"Task 1.4: Describe the characteristics of the study sample by salary status and predefine potential covariates for multiple regression model
table1(~ Rank + Prof + Discipline + Yrs.since.phd + Yrs.service + NPubs + Ncits + Salary + Sex | as.factor(high.salary), data = salary)| 0 (N=255) |
1 (N=142) |
Overall (N=397) |
|
|---|---|---|---|
| Rank | |||
| AssocProf | 63 (24.7%) | 1 (0.7%) | 64 (16.1%) |
| AsstProf | 67 (26.3%) | 0 (0%) | 67 (16.9%) |
| Prof | 125 (49.0%) | 141 (99.3%) | 266 (67.0%) |
| Prof | |||
| AssProfessor | 130 (51.0%) | 1 (0.7%) | 131 (33.0%) |
| Professor | 125 (49.0%) | 141 (99.3%) | 266 (67.0%) |
| Discipline | |||
| A | 127 (49.8%) | 54 (38.0%) | 181 (45.6%) |
| B | 128 (50.2%) | 88 (62.0%) | 216 (54.4%) |
| Yrs.since.phd | |||
| Mean (SD) | 19.1 (13.5) | 28.1 (9.29) | 22.3 (12.9) |
| Median [Min, Max] | 15.0 [1.00, 56.0] | 28.0 [11.0, 56.0] | 21.0 [1.00, 56.0] |
| Yrs.service | |||
| Mean (SD) | 14.9 (13.4) | 22.5 (10.7) | 17.6 (13.0) |
| Median [Min, Max] | 9.00 [0, 57.0] | 20.0 [2.00, 60.0] | 16.0 [0, 60.0] |
| NPubs | |||
| Mean (SD) | 17.6 (13.2) | 19.2 (15.2) | 18.2 (14.0) |
| Median [Min, Max] | 13.0 [1.00, 69.0] | 13.0 [1.00, 69.0] | 13.0 [1.00, 69.0] |
| Ncits | |||
| Mean (SD) | 38.6 (15.8) | 43.1 (18.6) | 40.2 (16.9) |
| Median [Min, Max] | 35.0 [1.00, 90.0] | 38.5 [1.00, 90.0] | 35.0 [1.00, 90.0] |
| Salary | |||
| Mean (SD) | 95000 (14600) | 147000 (20400) | 114000 (30300) |
| Median [Min, Max] | 96200 [57800, 120000] | 144000 [121000, 232000] | 107000 [57800, 232000] |
| Sex | |||
| Female | 30 (11.8%) | 9 (6.3%) | 39 (9.8%) |
| Male | 225 (88.2%) | 133 (93.7%) | 358 (90.2%) |
Task 1.5: Fit a multiple logistic regression to determine whether the chance of getting a high salary differed between male and female professors after adjusting for potential covariates. Interpret the finding.
Task 1.5.1: Simple Model
m.1 = glm(high.salary ~ Sex , family = "binomial", data = salary)
summary(m.1) ##
## Call:
## glm(formula = high.salary ~ Sex, family = "binomial", data = salary)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2040 0.3801 -3.168 0.00154 **
## SexMale 0.6782 0.3955 1.715 0.08636 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 517.75 on 396 degrees of freedom
## Residual deviance: 514.52 on 395 degrees of freedom
## AIC: 518.52
##
## Number of Fisher Scoring iterations: 4publish(m.1)## Variable Units OddsRatio CI.95 p-value
## Sex Female Ref
## Male 1.97 [0.91;4.28] 0.08636Task 1.5.2: Multiple Model
m.2 = glm(high.salary ~ Sex + Prof + Discipline + Yrs.since.phd + NPubs, family = "binomial", data = salary)
summary(m.2)##
## Call:
## glm(formula = high.salary ~ Sex + Prof + Discipline + Yrs.since.phd +
## NPubs, family = "binomial", data = salary)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.941036 1.140946 -5.207 1.92e-07 ***
## SexMale 0.215592 0.503442 0.428 0.668480
## ProfProfessor 5.031301 1.037844 4.848 1.25e-06 ***
## DisciplineB 0.997854 0.259625 3.843 0.000121 ***
## Yrs.since.phd 0.008645 0.012943 0.668 0.504147
## NPubs 0.004598 0.008735 0.526 0.598648
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 517.75 on 396 degrees of freedom
## Residual deviance: 363.49 on 391 degrees of freedom
## AIC: 375.49
##
## Number of Fisher Scoring iterations: 7publish(m.2)## Variable Units OddsRatio CI.95 p-value
## Sex Female Ref
## Male 1.24 [0.46;3.33] 0.6684800
## Prof AssProfessor Ref
## Professor 153.13 [20.03;1170.79] < 1e-04
## Discipline A Ref
## B 2.71 [1.63;4.51] 0.0001213
## Yrs.since.phd 1.01 [0.98;1.03] 0.5041467
## NPubs 1.00 [0.99;1.02] 0.5986477- Interpretation: There is no evidence (P= 0.64) that male professors were associated with different odds of getting a high salary than their female counterparts.
Task 2: Develop an algorithm to predict the chance of getting a high salary among professors.
Task 2.1: Study design
- The cross-sectional investigation of 397 professors was carried out to develop an algorithm for predicting the chance of getting a high salary.
Task 2.2: Import data set
Task 2.3: Describe characteristics of the study sample
Task 2.4: Check the distribution of professors’ salaries
Task 2.5: Check potential predictors of professors’ salaries
m.3 = bic.glm(high.salary ~ Sex + Prof + Discipline + Yrs.since.phd + Yrs.service + NPubs + Ncits, strict = F, OR = 20, glm.family = "binomial", data = salary)
summary(m.3)##
## Call:
## bic.glm.formula(f = high.salary ~ Sex + Prof + Discipline + Yrs.since.phd + Yrs.service + NPubs + Ncits, data = salary, glm.family = "binomial", strict = F, OR = 20)
##
##
## 9 models were selected
## Best 5 models (cumulative posterior probability = 0.9093 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 -5.978e+00 1.129453 -6.336e+00 -5.559e+00 -5.691e+00
## SexMale 4.6 1.053e-02 0.117232 . . .
## ProfProfessor 100.0 5.202e+00 1.018200 5.222e+00 5.198e+00 5.033e+00
## DisciplineB 100.0 9.733e-01 0.254090 9.797e-01 9.618e-01 1.006e+00
## Yrs.since.phd 5.3 4.565e-04 0.003544 . . 9.694e-03
## Yrs.service 2.1 1.670e-05 0.001626 . . .
## NPubs 4.7 1.025e-05 0.002167 . . .
## Ncits 51.8 9.540e-03 0.010718 1.837e-02 . .
##
## nVar 3 2 3
## BIC -1.993e+03 -1.993e+03 -1.988e+03
## post prob 0.425 0.408 0.027
## model 4 model 5
## Intercept -6.421e+00 -6.542e+00
## SexMale . 2.341e-01
## ProfProfessor 5.090e+00 5.206e+00
## DisciplineB 1.013e+00 9.837e-01
## Yrs.since.phd 7.581e-03 .
## Yrs.service . .
## NPubs . .
## Ncits 1.802e-02 1.843e-02
##
## nVar 4 4
## BIC -1.988e+03 -1.988e+03
## post prob 0.025 0.024
##
## 1 observations deleted due to missingness.imageplot.bma(m.3)
Task 2.6: Fit the model. Write a prediction algorithm.
m.4 = glm(high.salary ~ Prof + Discipline + Ncits, family = "binomial", data = salary)
summary(m.4)##
## Call:
## glm(formula = high.salary ~ Prof + Discipline + Ncits, family = "binomial",
## data = salary)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.335608 1.083334 -5.848 4.97e-09 ***
## ProfProfessor 5.221741 1.016653 5.136 2.80e-07 ***
## DisciplineB 0.979729 0.254899 3.844 0.000121 ***
## Ncits 0.018368 0.007614 2.412 0.015853 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 517.75 on 396 degrees of freedom
## Residual deviance: 358.41 on 393 degrees of freedom
## AIC: 366.41
##
## Number of Fisher Scoring iterations: 7Task 2.7: Develop a nomogram for implicating the prediction model
ddist = datadist(salary)
options(datadist = 'ddist')
m = lrm(high.salary ~ Prof + Discipline + Ncits, data = salary)
summary(m)## Effects Response : high.salary
##
## Factor Low High Diff. Effect S.E. Lower 0.95
## Ncits 28 50 22 0.4040900 0.16751 0.07577000
## Odds Ratio 28 50 22 1.4979000 NA 1.07870000
## Prof - AssProfessor:Professor 2 1 NA -5.2217000 1.01680 -7.21450000
## Odds Ratio 2 1 NA 0.0053979 NA 0.00073581
## Discipline - A:B 2 1 NA -0.9797300 0.25490 -1.47930000
## Odds Ratio 2 1 NA 0.3754100 NA 0.22779000
## Upper 0.95
## 0.732410
## 2.080100
## -3.228900
## 0.039599
## -0.480140
## 0.618700nom.m = nomogram(m, fun = function(x)1/(1+exp(-x)), fun.at = c(.001, .01, .10, seq(.2, .8, by = .2), .90, .99, .999), funlabel = "Predicted likelihood of a high salary")
plot(nom.m)