Bayesian Logistic Regression
Carolyn J. Anderson
11/14/2021
In this docuement at examples of both binary and multicategory (nominal) logistic regression. I will mostly rely on commments alreay include in an R script that I wrote.
In some cases I included mle code because I wasn’t sure if I was doing things correctly in brms.
The data set contains details about fatal police shooting from 2015 up through were downloaded from the Washington Post wesite. https://github.com/washingtonpost/data-police-shootings
Retrived 11/12/2021
The aim is to try to find out who is dying as a result of a police shooting, with attention given to the race of the victim. We’ll use race as the “outcome” or “response” variable and use some of the other variables to answer our question.
For reference, I pull this from https://www.census.gov and data are for July 19, 2019
| Classification | Percent population |
|---|---|
| Race and Hispanic Origin | |
| White alone, percent | 76.3% |
| Black or African American alone, percent(a) | 13.4% |
| American Indian and Alaska Native alone, percent(a) | 1.3% |
| Asian alone, percent(a) | 5.9% |
| Native Hawaiian and Other Pacific Islander alone, percent(a) | 0.2% |
| Two or More Races, percent | 2.8% |
| Hispanic or Latino, percent(b) | 18.5% |
| White alone, not Hispanic or Latino, percent | 60.1% |
Set up
Libraries
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, unionlibrary(tidyr)
library(lubridate) ##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(brms)## Loading required package: Rcpp## Loading 'brms' package (version 2.16.1). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:stats':
##
## arlibrary(coda)
library(nnet) Data
# Set up and look
setwd("D:/Dropbox/edps 590BAY/Lectures/11 Logistic Regression/Fatal Police Shootings Data")
df <- read.csv("fatal-police-shootings-data.csv", sep=",")
df <- as.data.frame(df)
#--- Recode date
df2 <- separate(df, date, c("month", "day", "year"))
#--- Only want pandemic year
# n=782, and 497 in November (so far)
df20 <- df2[which(df2$year==2020),]
#--- Race
table(df20$race)##
## A B H N O W
## 120 15 244 171 9 3 459# --- Only use B (black), H (hispanic) and W (white)
df <- df20[(df20$race=="H" | df20$race=="B" | df20$race=="W"),]
df$nonwhite <- ifelse(df$race=="W", 0,1)
#--- Age
table(df20$age)##
## 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 2 3 8 13 18 12 20 21 17 20 26 33 33 33 27 43 40 31 33 28 33 22 35 30 37 20
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
## 25 16 27 17 15 10 29 16 17 12 11 11 8 7 8 9 10 8 8 11 8 9 6 2 4 5
## 67 68 69 70 73 74 75 79 81 83 84 88
## 6 1 4 1 2 1 2 1 1 1 1 1# --- only non-children...drop 8 year old
df <- df[which(df$age>12),]
# --- and very elderly
df <- df[which(df$age<80),]Binary Logistic Regression
nonwhite ~ 1 + gender + age + body_camera
plot(df$age, jitter(df$nonwhite,.1), type="p", col="blue")
summary(glm(nonwhite ~ 1 + gender + age + body_camera,
data=df,
family=binomial))##
## Call:
## glm(formula = nonwhite ~ 1 + gender + age + body_camera, family = binomial,
## data = df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6936 -1.1107 -0.6377 1.1179 2.2340
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.251785 0.496500 0.507 0.61207
## genderM 1.177437 0.441274 2.668 0.00762 **
## age -0.042922 0.006317 -6.794 1.09e-11 ***
## body_cameraTRUE 0.548289 0.187140 2.930 0.00339 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1172.7 on 846 degrees of freedom
## Residual deviance: 1099.5 on 843 degrees of freedom
## AIC: 1107.5
##
## Number of Fisher Scoring iterations: 4options (mc.cores=parallel::detectCores ()) # Run on multiple cores
df$gender <- as.factor(df$gender)
df$body_camera <- as.factor(df$body_camera)
logreg1 <- brm(nonwhite ~ 1 + gender + age + body_camera,
data=df,
family = bernoulli(link = "logit"),
chains=4,
cores=4
) ## Compiling Stan program...## Start samplingplot(logreg1) 
# or
#logit1.mcmc <- as.mcmc(logreg1)
#traceplot(logit1.mcmc)
# etc
summary(logreg1)## Family: bernoulli
## Links: mu = logit
## Formula: nonwhite ~ 1 + gender + age + body_camera
## Data: df (Number of observations: 847)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.21 0.51 -0.81 1.19 1.00 4514 3139
## genderM 1.23 0.45 0.40 2.16 1.00 4376 3063
## age -0.04 0.01 -0.06 -0.03 1.00 4885 2775
## body_cameraTRUE 0.55 0.19 0.17 0.92 1.00 4320 3126
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(conditional_effects(logreg1), ask = FALSE)
df$agesq <- df$age*df$age
logreg2 <- brm(nonwhite ~ 1 + gender + age + agesq + body_camera,
data=df,
family = bernoulli(link = "logit"),
chains=4,
cores=4,
iter=3000
)
plot(logreg2)
# or
#logit2.mcmc <- as.mcmc(logreg1)
#traceplot(logit2.mcmc)
# etc
summary(logreg2)
plot(conditional_effects(logreg2), ask = FALSE)Multinomial
race ~ age
# --- MLE a a check ---
summary( mlogit1.mle <- multinom(race ~ age , data = df) ) ## # weights: 9 (4 variable)
## initial value 930.524609
## final value 833.215440
## converged## Call:
## multinom(formula = race ~ age, data = df)
##
## Coefficients:
## (Intercept) age
## H -0.3937137 0.00104456
## W -1.0413164 0.04516759
##
## Std. Errors:
## (Intercept) age
## H 0.3326162 0.00939056
## W 0.2747852 0.00738029
##
## Residual Deviance: 1666.431
## AIC: 1674.431# --- Bayesain ---
mlogit1 <- brm(race ~ age,
data=df,
family = "categorical",
iter=2000,
chains=4,
cores=4)## Compiling Stan program...## Start samplingplot(mlogit1)
summary(mlogit1)## Family: categorical
## Links: muH = logit; muW = logit
## Formula: race ~ age
## Data: df (Number of observations: 847)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## muH_Intercept -0.39 0.33 -1.06 0.24 1.00 3521 2670
## muW_Intercept -1.05 0.27 -1.59 -0.54 1.00 3414 2661
## muH_age 0.00 0.01 -0.02 0.02 1.00 3074 2745
## muW_age 0.05 0.01 0.03 0.06 1.00 3171 2621
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(mlogit1)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
race ~ 1 + gender + age + agesq + body_camera + signs_of_mental_illness
df$agesq <- df$age*df$age
summary(mlogit2.mle <- multinom(race ~ 1 + gender + age + agesq + body_camera + signs_of_mental_illness, data=df) )## # weights: 21 (12 variable)
## initial value 930.524609
## iter 10 value 819.137025
## final value 814.018586
## converged## Call:
## multinom(formula = race ~ 1 + gender + age + agesq + body_camera +
## signs_of_mental_illness, data = df)
##
## Coefficients:
## (Intercept) genderM age agesq body_cameraTRUE
## H -0.9570956 -1.241696 0.1072640 -0.001396861 -0.1144645
## W -0.7079400 -1.752675 0.1262161 -0.001072402 -0.6366606
## signs_of_mental_illnessTRUE
## H -0.09645451
## W 0.57462926
##
## Std. Errors:
## (Intercept) genderM age agesq body_cameraTRUE
## H 0.04593341 0.1720886 0.01590044 0.0002873503 0.2336668
## W 0.06866298 0.2162783 0.01613493 0.0002299720 0.2091575
## signs_of_mental_illnessTRUE
## H 0.1048466
## W 0.1483841
##
## Residual Deviance: 1628.037
## AIC: 1652.037mlogit2 <- brm(race ~ 1 + gender + age + agesq + body_camera + signs_of_mental_illness,
data=df,
family = "categorical",
iter=3000,
chains=4,
cores=4)## Compiling Stan program...## Start samplingplot(mlogit2)
summary(mlogit2, priors=TRUE, mc_se=TRUE, prob=.95)## Family: categorical
## Links: muH = logit; muW = logit
## Formula: race ~ 1 + gender + age + agesq + body_camera + signs_of_mental_illness
## Data: df (Number of observations: 847)
## Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
## total post-warmup draws = 6000
##
## Priors:
## Intercept_muH ~ student_t(3, 0, 2.5)
## Intercept_muW ~ student_t(3, 0, 2.5)
##
## Population-Level Effects:
## Estimate MCSE Est.Error l-95% CI u-95% CI Rhat
## muH_Intercept -0.85 0.02 1.34 -3.38 1.88 1.00
## muW_Intercept -0.45 0.02 1.14 -2.53 1.97 1.00
## muH_genderM -1.41 0.02 0.96 -3.44 0.35 1.00
## muH_age 0.11 0.00 0.05 0.01 0.21 1.00
## muH_agesq -0.00 0.00 0.00 -0.00 -0.00 1.00
## muH_body_cameraTRUE -0.12 0.00 0.24 -0.59 0.34 1.00
## muH_signs_of_mental_illnessTRUE -0.11 0.00 0.29 -0.70 0.48 1.00
## muW_genderM -2.00 0.02 0.86 -3.91 -0.59 1.00
## muW_age 0.13 0.00 0.04 0.05 0.20 1.00
## muW_agesq -0.00 0.00 0.00 -0.00 -0.00 1.00
## muW_body_cameraTRUE -0.64 0.00 0.21 -1.04 -0.21 1.00
## muW_signs_of_mental_illnessTRUE 0.58 0.00 0.23 0.15 1.05 1.00
## Bulk_ESS Tail_ESS
## muH_Intercept 3727 2864
## muW_Intercept 3598 2526
## muH_genderM 3171 2361
## muH_age 5192 4822
## muH_agesq 5259 4524
## muH_body_cameraTRUE 5152 4632
## muH_signs_of_mental_illnessTRUE 5072 4078
## muW_genderM 3282 2359
## muW_age 4739 4365
## muW_agesq 4705 4383
## muW_body_cameraTRUE 5114 4145
## muW_signs_of_mental_illnessTRUE 5378 4416
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(mlogit2)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Mulitlevel race ~ …. + (1 | state)
df$agesq <- df$age*df$age
mlogit3 <- brm(race ~ 1 + gender + age + agesq + body_camera + signs_of_mental_illness
+ (1 | state),
data=df,
family = "categorical",
iter=3000,
chains=4,
cores=4)## Compiling Stan program...## Start samplingplot(mlogit3)
summary(mlogit3, digits=5)## Family: categorical
## Links: muH = logit; muW = logit
## Formula: race ~ 1 + gender + age + agesq + body_camera + signs_of_mental_illness + (1 | state)
## Data: df (Number of observations: 847)
## Draws: 4 chains, each with iter = 3000; warmup = 1500; thin = 1;
## total post-warmup draws = 6000
##
## Group-Level Effects:
## ~state (Number of levels: 47)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(muH_Intercept) 1.44 0.29 0.95 2.09 1.00 2247 3553
## sd(muW_Intercept) 0.87 0.19 0.55 1.29 1.00 2175 3499
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## muH_Intercept -1.31 1.44 -4.15 1.58 1.00
## muW_Intercept 0.15 1.17 -2.06 2.53 1.00
## muH_genderM -1.43 0.95 -3.46 0.30 1.00
## muH_age 0.09 0.06 -0.02 0.20 1.00
## muH_agesq -0.00 0.00 -0.00 0.00 1.00
## muH_body_cameraTRUE -0.43 0.28 -0.97 0.12 1.00
## muH_signs_of_mental_illnessTRUE -0.17 0.31 -0.79 0.44 1.00
## muW_genderM -2.15 0.84 -3.96 -0.72 1.00
## muW_age 0.11 0.04 0.03 0.19 1.00
## muW_agesq -0.00 0.00 -0.00 0.00 1.00
## muW_body_cameraTRUE -0.79 0.24 -1.27 -0.32 1.00
## muW_signs_of_mental_illnessTRUE 0.66 0.23 0.21 1.11 1.00
## Bulk_ESS Tail_ESS
## muH_Intercept 4899 4500
## muW_Intercept 5104 4455
## muH_genderM 4858 3782
## muH_age 7154 4784
## muH_agesq 7207 4689
## muH_body_cameraTRUE 7919 4631
## muH_signs_of_mental_illnessTRUE 7659 5150
## muW_genderM 4823 3605
## muW_age 7317 4766
## muW_agesq 7205 4976
## muW_body_cameraTRUE 8877 5007
## muW_signs_of_mental_illnessTRUE 6763 5119
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pp_check(mlogit3)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Comparing the 3 multinomial models
# --- need these before can compare
mlogit1 <- add_criterion(mlogit1, criterion="loo")
mlogit1$criteria$loo##
## Computed from 4000 by 847 log-likelihood matrix
##
## Estimate SE
## elpd_loo -837.3 13.8
## p_loo 4.1 0.2
## looic 1674.6 27.5
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.mlogit2 <- add_criterion(mlogit2, criterion="loo")## Warning: Found 2 observations with a pareto_k > 0.7 in model 'mlogit2'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.mlogit2$criteria$loo ##
## Computed from 6000 by 847 log-likelihood matrix
##
## Estimate SE
## elpd_loo -826.8 15.2
## p_loo 12.8 1.3
## looic 1653.7 30.4
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 845 99.8% 3889
## (0.5, 0.7] (ok) 0 0.0% <NA>
## (0.7, 1] (bad) 2 0.2% 596
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.mlogit3 <- add_criterion(mlogit3, criterion="loo")
mlogit3$criteria$loo ##
## Computed from 6000 by 847 log-likelihood matrix
##
## Estimate SE
## elpd_loo -734.4 19.1
## p_loo 61.1 3.6
## looic 1468.8 38.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 839 99.1% 1359
## (0.5, 0.7] (ok) 8 0.9% 1810
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.loo(mlogit1, mlogit2, mlogit3, compare=TRUE, moment_match=TRUE, save_psis=TRUE)## Output of model 'mlogit1':
##
## Computed from 4000 by 847 log-likelihood matrix
##
## Estimate SE
## elpd_loo -837.3 13.8
## p_loo 4.1 0.2
## looic 1674.6 27.5
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'mlogit2':
##
## Computed from 6000 by 847 log-likelihood matrix
##
## Estimate SE
## elpd_loo -826.8 15.2
## p_loo 12.8 1.3
## looic 1653.7 30.4
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 845 99.8% 3889
## (0.5, 0.7] (ok) 0 0.0% <NA>
## (0.7, 1] (bad) 2 0.2% 596
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'mlogit3':
##
## Computed from 6000 by 847 log-likelihood matrix
##
## Estimate SE
## elpd_loo -734.4 19.1
## p_loo 61.1 3.6
## looic 1468.8 38.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 839 99.1% 1359
## (0.5, 0.7] (ok) 8 0.9% 1810
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## mlogit3 0.0 0.0
## mlogit2 -92.4 13.9
## mlogit1 -102.9 15.5# mlogit3 is the best
conditional_effects(mlogit3, categorical=T)
Some probability Plots from third Multinomial model
There is probably a better way to do this, but …
samples3<- as.data.frame(mlogit3, row.names = NULL, optional = TRUE)
means3 <- colMeans(samples3)
head(means3,n=13)## b_muH_Intercept b_muW_Intercept
## -1.313789442 0.152373857
## b_muH_genderM b_muH_age
## -1.428269461 0.091062558
## b_muH_agesq b_muH_body_cameraTRUE
## -0.001157535 -0.431278184
## b_muH_signs_of_mental_illnessTRUE b_muW_genderM
## -0.174595892 -2.147101807
## b_muW_age b_muW_agesq
## 0.111311557 -0.000835940
## b_muW_body_cameraTRUE b_muW_signs_of_mental_illnessTRUE
## -0.787158879 0.661446979
## sd_state__muH_Intercept
## 1.440107256#b_muH_Intercept b_muW_Intercept b_muH_genderM
#-1.3390132668 0.1304870133 -1.4095520449
#
# b_muH_age b_muH_agesq b_muH_body_cameraTRUE
#0.0916567049 -0.0011652295 -0.4249331808 #
#
#b_muH_signs_of_mental_illnessTRUE b_muW_genderM b_muW_age
#-0.1799139529 -2.1234545910 0.1109525408
#
#b_muW_agesq b_muW_body_cameraTRUE b_muW_signs_of_mental_illnessTRUE
#-0.0008288717 -0.7849902824 0.6588053271
par(mfrow=c(2,3))
# Use these to compute probabitlies...see if brms will do this
####### Male, Camera and Mental Illiness
agerange <- seq(from=15, to=80, length.out=150)
pwhite <- matrix(NA, nrow=150, ncol=2)
phispan <- matrix(NA, nrow=150, ncol=2)
pblack <- matrix(NA, nrow=150, ncol=2)
i <- 1
for (age in agerange) {
pwhite[i,1] <- age
linpred1 <- means3[2] + means3[8] + means3[9]*age + means3[10]*age*age + means3[11] + means3[12]
phispan[i,1] <- age
linpred2 <- means3[1] + means3[3] + means3[4]*age + means3[5]*age*age + means3[6] + means3[7]
pblack[i,1] <- age
linpred3 <- 0
denmoinator <- exp(linpred1) + exp(linpred2) + exp(linpred3)
pwhite[i,2] <- exp(linpred1)/denmoinator
phispan[i,2] <- exp(linpred2)/denmoinator
pblack[i,2] <- exp(linpred3)/denmoinator
i <- i + 1
}
plot(pblack[,1], pblack[,2],type="l",col="black", lwd=2,
ylim=c(0,1),
xlab="Age",
ylab="Probability Victim",
main="Male, Camera and Mental Illiness")
lines(pwhite[,1], pwhite[,2],col="cyan", lwd=2)
lines(phispan[,1], phispan[,2],col="red", lwd=2)
##### Male, No Camera and Mental Illiness"
agerange <- seq(from=15, to=80, length.out=150)
pwhite <- matrix(NA, nrow=150, ncol=2)
phispan <- matrix(NA, nrow=150, ncol=2)
pblack <- matrix(NA, nrow=150, ncol=2)
i <- 1
for (age in agerange) {
pwhite[i,1] <- age
linpred1 <- means3[2] + means3[8] + means3[9]*age + means3[10]*age*age + means3[12]
phispan[i,1] <- age
linpred2 <- means3[1] + means3[3] + means3[4]*age + means3[5]*age*age + means3[7]
pblack[i,1] <- age
linpred3 <- 0
denmoinator <- exp(linpred1) + exp(linpred2) + exp(linpred3)
pwhite[i,2] <- exp(linpred1)/denmoinator
phispan[i,2] <- exp(linpred2)/denmoinator
pblack[i,2] <- exp(linpred3)/denmoinator
i <- i + 1
}
plot(pblack[,1], pblack[,2],type="l",col="black", lwd=2,
ylim=c(0,1),
xlab="Age",
ylab="Probability Victim",
main="Male, No Camera and Mental Illiness")
lines(pwhite[,1], pwhite[,2],col="cyan", lwd=2)
lines(phispan[,1], phispan[,2],col="red", lwd=2)
##### Male, No Camera and No Mental Illiness"
agerange <- seq(from=15, to=80, length.out=150)
pwhite <- matrix(NA, nrow=150, ncol=2)
phispan <- matrix(NA, nrow=150, ncol=2)
pblack <- matrix(NA, nrow=150, ncol=2)
i <- 1
for (age in agerange) {
pwhite[i,1] <- age
linpred1 <- means3[2] + means3[8] + means3[9]*age + means3[10]*age*age
phispan[i,1] <- age
linpred2 <- means3[1] + means3[3] + means3[4]*age + means3[5]*age*age
pblack[i,1] <- age
linpred3 <- 0
denmoinator <- exp(linpred1) + exp(linpred2) + exp(linpred3)
pwhite[i,2] <- exp(linpred1)/denmoinator
phispan[i,2] <- exp(linpred2)/denmoinator
pblack[i,2] <- exp(linpred3)/denmoinator
i <- i + 1
}
plot(pblack[,1], pblack[,2],type="l",col="black", lwd=2,
ylim=c(0,1),
xlab="Age",
ylab="Probability Victim",
main="Male, No Camera and No Mental Illiness")
lines(pwhite[,1], pwhite[,2],col="cyan", lwd=2)
lines(phispan[,1], phispan[,2],col="red", lwd=2)
####### Woman, Camera and Mental Illiness
agerange <- seq(from=15, to=80, length.out=150)
pwhite <- matrix(NA, nrow=150, ncol=2)
phispan <- matrix(NA, nrow=150, ncol=2)
pblack <- matrix(NA, nrow=150, ncol=2)
i <- 1
for (age in agerange) {
pwhite[i,1] <- age
linpred1 <- means3[2] + means3[9]*age + means3[10]*age*age + means3[11] + means3[12]
phispan[i,1] <- age
linpred2 <- means3[1] + means3[4]*age + means3[5]*age*age + means3[6] + means3[7]
pblack[i,1] <- age
linpred3 <- 0
denmoinator <- exp(linpred1) + exp(linpred2) + exp(linpred3)
pwhite[i,2] <- exp(linpred1)/denmoinator
phispan[i,2] <- exp(linpred2)/denmoinator
pblack[i,2] <- exp(linpred3)/denmoinator
i <- i + 1
}
plot(pblack[,1], pblack[,2],type="l",col="black", lwd=2,
ylim=c(0,1),
xlab="Age",
ylab="Probability Victim",
main="Woman, Camera and Mental Illiness")
lines(pwhite[,1], pwhite[,2],col="cyan", lwd=2)
lines(phispan[,1], phispan[,2],col="red", lwd=2)
####### Woman, no Camera and Mental Illiness
agerange <- seq(from=15, to=80, length.out=150)
pwhite <- matrix(NA, nrow=150, ncol=2)
phispan <- matrix(NA, nrow=150, ncol=2)
pblack <- matrix(NA, nrow=150, ncol=2)
i <- 1
for (age in agerange) {
pwhite[i,1] <- age
linpred1 <- means3[2] + means3[9]*age + means3[10]*age*age + means3[12]
phispan[i,1] <- age
linpred2 <- means3[1] + means3[4]*age + means3[5]*age*age + means3[7]
pblack[i,1] <- age
linpred3 <- 0
denmoinator <- exp(linpred1) + exp(linpred2) + exp(linpred3)
pwhite[i,2] <- exp(linpred1)/denmoinator
phispan[i,2] <- exp(linpred2)/denmoinator
pblack[i,2] <- exp(linpred3)/denmoinator
i <- i + 1
}
plot(pblack[,1], pblack[,2],type="l",col="black", lwd=2,
ylim=c(0,1),
xlab="Age",
ylab="Probability Victim",
main="Woman, no Camera and Mental Illiness")
lines(pwhite[,1], pwhite[,2],col="cyan", lwd=2)
lines(phispan[,1], phispan[,2],col="red", lwd=2)
####### Woman, no Camera and no Mental Illiness
agerange <- seq(from=15, to=80, length.out=150)
pwhite <- matrix(NA, nrow=150, ncol=2)
phispan <- matrix(NA, nrow=150, ncol=2)
pblack <- matrix(NA, nrow=150, ncol=2)
i <- 1
for (age in agerange) {
pwhite[i,1] <- age
linpred1 <- means3[2] + means3[9]*age + means3[10]*age*age
phispan[i,1] <- age
linpred2 <- means3[1] + means3[4]*age + means3[5]*age*age
pblack[i,1] <- age
linpred3 <- 0
denmoinator <- exp(linpred1) + exp(linpred2) + exp(linpred3)
pwhite[i,2] <- exp(linpred1)/denmoinator
phispan[i,2] <- exp(linpred2)/denmoinator
pblack[i,2] <- exp(linpred3)/denmoinator
i <- i + 1
}
plot(pblack[,1], pblack[,2],type="l",col="black", lwd=2,
ylim=c(0,1),
xlab="Age",
ylab="Probability Victim",
main="Woman, no Camera and no Mental Illiness")
lines(pwhite[,1], pwhite[,2],col="cyan", lwd=2)
lines(phispan[,1], phispan[,2],col="red", lwd=2)