Clayton Besaw is a PhD student in the Political Science department at the Univeristy of Central Florida. This site details his research and thoughts concerning politics, data-science, and everything in-between.
This is a repost of something that I shared on social media. Since this blog is about politics and data, I thought that it was appropriate.
During these times of magnified focus on extremist acts, we need to keep in mind some ideas about extremist and terrorist violence.
1) Overall, terrorism is a rare phenomenon and the statistical likelihood of you being killed or attacked is small.^{1}
2) The goal of terrorist acts is not to kill you, but to make you angry or fearful. They want you to overreact or give up. These are their strategies and it is the same for Right-wing terrorists, Left-wing terrorists, or Islamic terrorists.^{2}
3) The majority of those that commit such acts of violence are not psychopathic or hold un-malleable beliefs of hatred. Extremism is a social phenomenon that is the product of many socio-psychological, economic, political, and social factors.^{3}
We should be most cognizant of this third point. Extremist groups, whether they are Islamist or Nazi, are often like gangs. They produce intense peer pressure on group members to join and to maintain their status as a member. Their control over members is strengthened by extreme othering of outsiders. This othering thrives on commonly held arguments that paint in-group members as being psychologically or morally defunct.^{4}
So what can we do? Sweden and Germany have had success with creating accessible and understanding decradicalization programs. One such program is called EXIT.^{5} EXIT Deutschland has been instrumental in helping members of violent Nazi gangs and radical Islamist groups quit their lives as members of these organizations. EXIT sets up anonymous hot lines that provide advice, counseling, and understanding. If a member wishes to quit their group, EXIT can help to provide former extremists with safe housing and protection from reprisal.
Americans from the left, right, and center have fear of those that have taken up the mantel of extremism, but fear is our own worst enemy. We owe it to ourselves to empathize and understand those that commit such acts and to provide a bridge for them to escape the siren call of extremism and the illusion of safety and comradeship that it creates.
As you think of ways to make a difference politically or socially (whatever that may be), we can all agree that our democratic duty is to help those that are trapped in cycles of violence and extremist social pressure. Here are some organizations that one may think of donating to or volunteering with in the United States. Most importantly as well. If you personally know someone who you fear is courted or already a part of an extremist organization take time to reach out to them and to point them in the direction of these organizations. Even if they do not respond, you have planted the seed of doubt. If their doubt grows, these organizations can help them live again.
Kydd, Andrew H., and Barbara F. Walter. “The strategies of terrorism.” International Security 31, no. 1 (2006): 49-80. ↩
Victoroff, Jeff. “The mind of the terrorist A review and critique of psychological approaches.” Journal of Conflict resolution 49, no. 1 (2005): 3-42. ↩
Atran, Scott, Hammad Sheikh, and Angel Gomez. “Devoted actors sacrifice for close comrades and sacred cause.” Proceedings of the National Academy of Sciences 111, no. 50 (2014): 17702-17703. ↩
Lombardi, Marco, Eman Ragab, and Vivienne Chin, eds. Countering radicalisation and violent extremism among youth to prevent terrorism. Vol. 118. IOS Press, 2014. ↩
This post chronicles my attempt at the Kaggle Titanic competition. The objective of this competition is to build a classifier that correctly predicts passenger survival in the training data. My high score as a result of this effort is 370th out of ~5800, or roughly the top 7%. This post is only part 1 of this process. It outlines the use of logistic regression as a classifier and model specification. Part 2 will detail the random forest approach, which had much better results.
Data Munging and Variable Inspection
First, import the CSV data anyway you like.
train <- read.csv(".../train.csv")
test <- read.csv(".../test.csv")
For this process as well, let’s combine the training and test sets for easier data manipulation. The test set starts at passenger 892, so we can resplit the data later for analysis.
test$Survived <- NA #test needs matching column length for binding
comb_df <- rbind(train, test)
Use the head() function to inspect the first six observations for each variable. Also, the str() function is useful for looking at data-frame component characteristics. These include data type, measurement levels, and matching variable names for this information.
head(comb_df)
## PassengerId Survived Pclass
## 1 1 0 3
## 2 2 1 1
## 3 3 1 3
## 4 4 1 1
## 5 5 0 3
## 6 6 0 3
## Name Sex Age SibSp
## 1 Braund, Mr. Owen Harris male 22 1
## 2 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female 38 1
## 3 Heikkinen, Miss. Laina female 26 0
## 4 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35 1
## 5 Allen, Mr. William Henry male 35 0
## 6 Moran, Mr. James male NA 0
## Parch Ticket Fare Cabin Embarked
## 1 0 A/5 21171 7.2500 S
## 2 0 PC 17599 71.2833 C85 C
## 3 0 STON/O2. 3101282 7.9250 S
## 4 0 113803 53.1000 C123 S
## 5 0 373450 8.0500 S
## 6 0 330877 8.4583 Q
One should quickly notice that the the dependent variable (survived) is not a factor. Pclass also appears to be misclassified as an integer as well. These should be transformed into the factor data type with ordered levels. I make new variables for each factor because I prefer leaving the original variables in their original form.
Now let’s examine the extent of missing data. This can be easily done by using the VIM package and the aggr plot function. Let’s also quickly remove the survived variable since it is obviously missing for the test data.
##
## Variables sorted by number of missings:
## Variable Count
## Age 0.2009167303
## Fare 0.0007639419
## PassengerId 0.0000000000
## Pclass 0.0000000000
## Name 0.0000000000
## Sex 0.0000000000
## SibSp 0.0000000000
## Parch 0.0000000000
## Ticket 0.0000000000
## Cabin 0.0000000000
## Embarked 0.0000000000
## Pclass2 0.0000000000
The histogram and resulting table indicate that roughly 20% of age values are missing, and that a very small percent (0.0007) of the fare values are missing.^{1} Logistic regression can handle missing values through listwise deletion, but this is not a sound strategy if we have the tools to impute missing values. Methods like Random Forests can’t work with missing values anyways, so we should really find a nice imputation solution. Let’s assume that these values are missing at random and are of an adequately low percentage of full data-frame representation.^{2}
Since fare is only missing one, it is easily enough to impute using some central tendency measure. The mean is 33.29, while the median is 14.45. Inspecting the fare variable shows a number of problems.
summary(comb_df$Fare)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 7.896 14.450 33.300 31.280 512.300 1
First, there are Fares of 0. Did these people get on for free? Probably not. Second, it is clear upon inspection that outliers are pushing the mean Fare value up. Now we have more than just a simple imputation to deal with. We can impute the missing Fare with the median fare based on his class using the Rpart package, but first let’s set the zero values to NA and impute median values based on class.
comb_df$fare2 <- ifelse(comb_df$Fare==0, NA, comb_df$Fare)
summary(comb_df$fare2) #inspect to see no more zero value, and now 18 NAs
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 3.171 7.925 14.500 33.730 31.330 512.300 18
library(rpart)
fareimp <- rpart(fare2 ~ Pclass, data=comb_df[!is.na(comb_df$fare2),],
method="anova")
comb_df$fare2[is.na(comb_df$fare2)] <- predict(fareimp, comb_df[is.na(comb_df$fare2),])
summary(comb_df$fare2) #inspect to see no more NAs
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.171 7.925 14.500 33.880 31.390 512.300
Now let’s impute age based on multiple covariates.
ageimp <- rpart(Age ~ Pclass2 + Sex + SibSp + Parch + fare2 + Embarked,
data=comb_df[!is.na(comb_df$Age),],
method="anova")
comb_df$age2 <- NA
comb_df$age2[is.na(comb_df$age2)] <- predict(ageimp, comb_df[is.na(comb_df$age2),])
summary(comb_df$age2) #inspect to make sure.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 9.60 27.43 27.43 29.63 32.10 40.55
Finally in terms of data extraction, let’s look at titles housed in the name variable. People more clever than I noticed that titles were consistent throughout the names listed in the passenger data. Historically, titles were also associated with different personal characteristics such as gender, class, nobility, etc.
comb_df$Name <- as.character(comb_df$Name)
strsplit(comb_df$Name[2], split='[,.]') #want to isolate the "title" string segment
Finally, let’s make a women and children first variable. This is based on the common anecdotal story of life boats only accepting women and children before adult males.
We have done a lot of work here! Let’s run some models. First resplit the combined data set back into training and test sets.
train <- comb_df[1:891,]
test <- comb_df[892:1309,]
Logistic Regression
Logistic regression is a classic classifier in statistics and econometrics. It is also heavily used in political science and sociology as a workhorse modeling method. Like linear regression and classification methods, we are interested in understanding the probability of an outcome \(Y=1\) occurring conditional on a set of \(n\) predictor variables within \(x\prime\).
\[ Pr(Y = 1|x\prime) \]
One can read the above as being the probability of \(Y\) equaling 1 given the measures found in \(x\prime\). To model this predictive relationship, we can utilize conventions from linear regression.
First we can represent the above conditional probability as:
\[ Pr(Y = 1|x\prime) = p(x\prime) \]
Next, propose a linear function of this probability:
Unfortunately, stopping here would not help us. Our model would sometimes result in probabilities that are \(p(x\prime) < 0\) or \(p(x\prime) > 1\). It is impossible to have probabilities outside of the range \(p \in [0,1]\). Fortunately, we can use the logistic function to guarantee outputs between 0 and 1 as such:
Reiterating basic probability, the odds tells us about the probability of \(Y = 1\). Odds close to 0 indicate very low probabilities for \(Y = 1\), while probabilities ever approaching \(\infty\) indicate very high probabilities for \(Y = 1\). By logging both sides of the above odds equation, we can calculate the logit or logged odds:
This is the basis for using logistic regression to make classification predictions. We can interpret the influence of predictor \(x_{i}\) as the change in the logged odds of \(Y = 1\) by \(\beta_{i}\). This is not a linear increase like the linear regression interpretation, but is a conditional increase based on the current value of \(x_{i}\). Even though the relationship is non-linear, we can generalize the interpretation of \(\beta_{i}\). If \(\beta_{i}\) is positive, then increasing the value of \(x_{i}\) will increase the logged odds of \(Y = 1\). As such, a negative \(\beta_{i}\) would have the opposite effect.
The above seems easy enough, but how can we know the values of \(\beta_{0}\), \(\beta_{i}\), etc..? These parameters are generally estimated though the maximum likelihood method. Let’s assume \(n\) observations and that \(x\) and \(\beta\) are sets of covariates and their corresponding beta estimates. First start with the likelihood function:
Finally, one can find the maximum likelihood estimates for a specific parameter by taking the partial derivative of the log likelihood function. Say we wanted the estimate \(\beta_{j}\) from the set of estimates \(\beta\). We take the partial derviative of the log likelihood function with respect to \(\beta_{j}\)
From here, we have the basis for estimating \(\beta_{j}\). The resulting beta values are estimates because the partial derivative function is transcendental.
Logistic Regression using GLM
So let’s put the above into action. Logitistic regression models are easy to implement in R using the glm() function. Here I use the logit link.
Multiple imputation assumes that missing values are missing at random (MAR). One should be careful when performing multiple imputation on real world data, as the MAR assumption is often not feasible. I am not sure why Age has so many missing values here. It could be a design by Kaggle, which is what I personally assume for this analysis. Normally, if one variable has such a high degree of missingness, it is unlikely to be missing at random. In addition, missing values should normally not be more than 5%-10% of a measure’s representation. Having much higher degrees of missingness is often a sign of MAR violation. Consequently, with more missingness, you are more likely to artificially construct patterns within the imputed data if your baseline data is unrepresentative of the real population.↩
Here is some python code for making a simple number guessing game.
importrandomguessesTaken=0print("Hello! What is your name?")my_name=input("Name: ",)number=random.randint(1,20)print("Well, "+my_name+", I am thinking of a number between 1 and 20.")whileguessesTaken<6:print("Take a guess.")guess=int(input())guessesTaken=guessesTaken+1ifguess<number:print("Your guess is too low!")ifguess>number:print("Your guess is too high!")ifguess==number:breakifguess==number:guessesTaken=str(guessesTaken)print("Good job, "+my_name+"! You guessed my number in "+guessesTaken+' guesses!')ifguess!=number:number=str(number)print("Nope. The number I was thinking of was "+number)
This post is really only a test for Rmarkdown integration into Jekyll. Substantively, it is a first lab lesson for applied Bayesian estimation in R.
Load the LearnBayes^{1} and Lattice packages.
1. Please replicate the sleeping habits analysis using a discrete prior (in the session 1-2)
Setting up a discrete prior
This code is used to create a a discrete prior concerning our beliefs about the proportion college student sleeping habits (p):
The code above is two-fold, first, it sets up the range of possible proportion values for the amount of students that could be sleeping at least eight hours a night. Second, it creates our discrete prior by specifying our weighted beliefs about each distinct proportion (p).
Next, we plot the our priors for each proportion (p) by using the lot function:
The histogram above visualizes our prior belief in each value of p. The original guess of .3 (from lesson 1-2) being the most probable is highlighted clearly in the histogram. We see that our prior beliefs are greatest between .2 and .4, with.3 being the peak in the middle.
Evidence and Posterior
Next, we want to calculate the posterior probabilities for p by utilizing the data on college student sleeping habits. In the example study, only 11 out of 27 participants were found to have slept a sufficient number of hours per night. By using the pdisc function from the LearnBayes package, we can utilize our established priors to compute our posterior probabilities. The code below is used to create our posterior probabilities:
First, the code above establishes our evidence (data) from the example study (11 students out of 27). Second, the posterior values are calculated using the pdisc function, which takes our discrete values (p), our priors, and the evidence. The function below provides simple table of the work conducted so far:
Here we can easily see the prior and posterior probabilies for each level of p.
Comparsion of Prior and Posterior
By using the lattice package, and the xyplot function, R can visualize the comparison between the prior and posterior probabilities. The code below sets up the comparison graphs:
The above code creates a dataset composed of the prior and posterior distributions. Next, this dataset is used to create a set of comparison graphs for the priors and the posteriors:
Here we can see a direct visual comparison between the prior and posterior probability distributions. Taking the evidence into account, we see that posterior probability peak shifts closer to .4. Overall, the main set of probability distributions is settled in-between .25 and .45 in terms of proportion (p).
2. Please replicate the sleeping habits analysis using a beta prior (in the session 1-2b).
Setting up the Beta Prior
The specification of a beta distribution is useful for bayesian analysis when it concerns the behavior of percentages and proportions.^{2} The beta distribution is reflected by two shape parameters a and b, which can represent our prior beliefs about the proportion of students who get at least eight hours of sleep a night p. The code below creates the shape parameters for the beta prior by using the beta.select function, which utilizes our examples prior belief about two quantiles of p, the median (.3) and and 90th percentile (.5).^{3} Our example study believes that the median is equally likely (.5) to be smaller or larger than the median (.3), and that p is less than .5 (90th Percentile), with 90% confidence.
The beta.select function provides us with the shape parameters of a(3.26) and b(7.19). This was done by first creating two variables that lists our beliefs about the median (.3) and 90th percentile (.5) quantiles, and using them as arguments in the beta.select function.
Evidence and Posterior
For the evidence and likelihood function, we retain the previous example study data of 11 college students out of 27 who receive at least eight hours of sleep. A conjugate analysis is possible with the beta prior, because once you combine the beta prior with the evidence, a beta posterior is computed. The code below sets up the beta posterior by using the dbeta command:
Here we can see a visualized comparison between the beta prior and posterior. As with the discrete posterior, the peak probability density for the beta posterior shifts towards .4 for p, after we take the study evidence into account.
Analyzing the beta posterior
This section looks at three different ways in which the beta posterior can be analyzed. For this, we ask the question; Is it likely that the proportion of heavy sleepers is greater than .5?
1) The first way to analyze this question is to compute the probability of p being greater than .5 by using the pbeta function, which is the beta distribution function. The code below calculates the probability of our question being right:
The probability calculated is 6.9%, which is small. We can be confident that p is probably not more than .5, or more than half of college student sleepers.
2) The seoncd way to analyze this question is to create a credible interval estimate, which is a bayesian take on frequentist confidence intervals.^{4} Credible Intervals can be calculated by using the qbeta function in R. The code below calculates the credible intervals for the beta posterior:
The calculated 90% credible interval for our beta posterior is (0.256, 0.513). Again, it seems unlikely that p is greater than .5 because the credible interval only reaches past .5 on the uppertail, with the majority of the interval being below .5.
3) The third way to analyze this question is to use a simulation of a large number of values from the beta posterior. To do this, the rbeta function is used to generate a large amount of random proportion values based on the beta shape paramters. The code below computes the simulation and visualizes the results using the hist function:
The histogram above shows newly simulated sample for p. Next we can replicate the previous two posterior analyzes by using the sum and quantile functions. The code below calculates the probability and credible interval for the question of interest:
Both the probability (6.9% vis a vis 6.7%) and the credible interval estimates (0.25, 0.51 vis a vis .26, .52) are very close to the original exact values of the beta posterior. The results still lead to the analyst to conclude that p is unlikely to be greater than .5, but what if we can increase percision by greater simulation? The next block of code replicates the previous simulation procedure but with 1,000,000 generated values for p instead.
After increasing the simulation to a 1,000,000 values of p, we find even closer estimates of the probability of p > .5 (6.9% vis a vis 6.8%) and the credible interval (0.25, 0.51 vis a vis .25, 0.51). We can be fairly confident in our previous findings concerning the question. After taking into account the study evidence, it is unlikely that the proportion of college students who have at least eight hours of sleep a night is greater than .5.
3. Please replicate the birth rate comparison analysis using a Poisson Likelihood (in the session 1-5).
For the study of birth rates amongst women with, and without college degrees, we utilize the Poisson model. Birthrates qualify as count data, which can be better analyzed by using poisson and negative binomial modeling.
Setting up the Prior
For a Poisson likelihood model, one can use a gamma conjugate prior. Again, this type of prior only needs to model our beliefs about two parameters, a(occurances) and b(intervals). The code below sets up our gamma prior.
Evidence and Posterior
Next, we wish to take into account the evidence about birthrates amongst women without college degrees, and with college degrees. This data comes from the General Social Survey (GSS), which interviewed women on educational attainment and number of children. From here, group 1 defines women without a college degree, and group 2 defines women with a college degree. The code below establishes the evidence (data) for both groups of women, and then finds the 95% credible interval for number of children for each group.
The posterior mean number of children born to group 1 is 1.95, the mode is 1.94, and the 95% credible interval is (1.70, 2.22). The code below computes the posterior data for group 2.
The posterior mean number of children born to group 2 is 1.51, the mode is 1.48, and the 95% credible interval is (1.17, 1.89).
Posterior Group Comparison
Comparing the posterior results between group 1 (without degrees) and group 2 (with degrees) shows that women with degrees appear to have less children (y) on average than their counterparts without degrees (1.51 versus 1.95). The posterior credible intervals also show two mostly different ranges for number of children (y) had by both groups of women (1.70, 2.22 versus 1.17, 1.89). The code below creates a visualized comparison of the gamma prior, posterior for group 1, and posterior for group 2:
The plot above confirms the previous observation about the posterior comparison between the two groups of women. Women without degrees peaks around 2 children, while women with degrees peaks at around 1.5 children.
Birth Rate Prediction
To make predictions about the number of children women will have depending on their college degree status, we can utilize the posteriors found above with a negative binomial distribution. This can be done in r by using the dnbinom function. The code below creates predictions of birth rates (y) for both groups of women:
The code above creates predicted probabilities for number of children born on an interval of 0 to 10. The dnbinom function uses two arguments to create the predicted probabilities. The first is the size, or dispersion paramters, and the second argument is the mean.
The code below creates a visualized comparison of the predicted probabilities for both groups of women:
As expected, the plot above shows that women without degrees have higher predicted probabilities to have two or more children. In turn, women with a degree have relatively higher predicted probabilities for having a single child or no child, when compared to women without a degree.
LearnBayes is a package created by Jim Albert at Bowling Green University. Documentation can be found at http://cran.r-project.org/web/packages/LearnBayes/LearnBayes.pdf. ↩
Kruschke, John. Doing Bayesian data analysis: A tutorial introduction with R. Academic Press, 2010. ↩
Beta priors can also be established thorugh the use of our beliefs concerning the beta distribution’s mean and standard deviation, but this can be problematic, because the analyst must make sure the standard deviation makes sense for a beta distribution (Kruschke, J, 2010). ↩
Credible Intervals are essentially analagous to Confidence Intervals in frequentist statistical analysis, but differ in a key way. Credible intervals incorporate contextual information from the prior, while confidence intervals take only the data (evidence) into account. For more on credible intervals vis a vis confidence intervals, please see: Greenberg, Edward. Introduction to Bayesian econometrics. Cambridge University Press, 2012. ↩
I have finally figured out Github page hosting with Jekyll. I would like to thank Barry Clark for his incredible walkthrough for using Jekyll Now as a static site generator. Additionally, I would like to acknowledge Mark Otto as the author/creator of the Lanyon theme for Jekyll that is used by this site.