While scientists estimate that about 44,000 kilograms of meteoric material falls on Earth each day only a tiny fraction is recovered and classified. I find these landings provide valuable information across different fields, offering important scientific and historical insights that help us understand more about our universe and our own history. One of the most important factors when studying meteorites is knowing whether they were a “fall” or a “find.” Determining this helps us assess whether the meteorite’s composition is still intact or if it has been deteriorated by environmental factors such as wind, water heat or biological factors like plants, bacteria, and animals. The elaboration of this decision tree not only provides information regarding the discovery status and year but it also accounts for the meteorite’s physical properties such as size. Knowing this status is critical as it allows us to extract valuable information regarding the meteorite’s true composition and history. Link to the data here

glimpse(meteorite)
## Rows: 45,716
## Columns: 10
## $ name        <chr> "Aachen", "Aarhus", "Abee", "Acapulco", "Achiras", "Adhi K…
## $ id          <int> 1, 2, 6, 10, 370, 379, 390, 392, 398, 417, 423, 424, 425, …
## $ nametype    <chr> "Valid", "Valid", "Valid", "Valid", "Valid", "Valid", "Val…
## $ recclass    <chr> "L5", "H6", "EH4", "Acapulcoite", "L6", "EH4", "LL3-6", "H…
## $ mass..g.    <dbl> 21, 720, 107000, 1914, 780, 4239, 910, 30000, 1620, 1440, …
## $ fall        <chr> "Fell", "Fell", "Fell", "Fell", "Fell", "Fell", "Fell", "F…
## $ year        <int> 1880, 1951, 1952, 1976, 1902, 1919, 1949, 1814, 1930, 1920…
## $ reclat      <dbl> 50.77500, 56.18333, 54.21667, 16.88333, -33.16667, 32.1000…
## $ reclong     <dbl> 6.08333, 10.23333, -113.00000, -99.90000, -64.95000, 71.80…
## $ GeoLocation <chr> "(50.775, 6.08333)", "(56.18333, 10.23333)", "(54.21667, -…

Describing the dataset

The original dataset comprised 45,000 observations however, due to a severe class imbalance, 97% “Found” versus only 3% “Fell”, it was necessary to extract a more representative sample for model training.

table(meteorite$fall)
## 
##  Fell Found 
##  1107 44609

To solve this I created a more balanced training set of 1,004 observations, consisting of 502 “Fell” and 502 “Found” records.

table(balanced_data$fall)
## 
##  Fell Found 
##   502   502

I also reorganized the data to improve clarity and model performance, years were grouped into centuries, latitudes and longitudes were converted into geographic regions, and mass was transformed using a logarithmic scale to address the extreme range in values. These adjustments resulted in a more intuitive and readable dataset that keeps the original information while significantly improving the model’s ability to learn from the data.

balanced_data %>%
  select(5,(last_col() - 3):last_col()) %>%
  glimpse()
## Rows: 1,004
## Columns: 5
## $ fall                <chr> "Fell", "Fell", "Fell", "Fell", "Fell", "Fell", "F…
## $ Century_Num         <dbl> 19, 21, 20, 21, 20, 20, 20, 20, 19, 19, 20, 20, 19…
## $ Log_Mass            <dbl> 3.000, 2.225, 2.903, 4.230, 1.884, 3.279, 4.394, 5…
## $ Period              <chr> "Industrial Revolution Era", "Digital Era", "WW I/…
## $ Geographical_Region <fct> Europe, Africa, Asia, Africa, Europe, Asia, Ocean/…

Splitting the dataset

Before fitting the tree, we will randomly split the dataset into two: train data and test data. The train data is used to learn the model and the test data is used to assess the prediction power of the model. The common practice is to split the data 80/20, 80 percent of the data serves to train the model, and 20 percent to make predictions. As the splitting is random, we will fix the seed number in R to ensure we all get the same results

set.seed(2026)
n = length(balanced_data$fall) # count number of rows to know the sample size
nt = round(n*0.80)  # compute the 80% of the sample size and round to an integer number
index = sample(1:n, nt) # randomly select "nt" numbers from 1 to n and storage them in the vector "index"
train=balanced_data[index,]  # take only the rows = index for the training data
test=balanced_data[-index,] # take all the rows that are not in "index" for the testing data

Fitting the decision tree

Once the data has been split into two, we will train the decision tree model using the rpart() function on the train dataset. It is necessary to specify what is the binary dependent variable, and what are the predictors. Additionally, the data source must be defined and the method argument set to class to ensure proper handling for this classification problem

mytree <- rpart(
  fall ~ Period+Century_Num+Log_Mass+year+Geographical_Region,  
  data = train, 
  method = "class"
)

Once the model has been fitted we can visualize it by plotting the decision tree:

rpart.plot(mytree)

Interpreting the nodes

Starting from the root node (the one at the top) the top label indicates the most frequent classification in that group (Fell or Found), the second number is the probability of the Found class and the third number is the percentage of the total observations contained in that node. From here the tree splits at year < 1973. The left path is for meteorites recorded before 1973 corresponding to 42% of the total observations. Among these, the probability of being Found is only 7%, leaving the 93% left as the Fell dominant class. Following the right path for meteorites recorded in 1973 or later, here the probability of being Found is 79% making it the predicted class. Analyzing the subsequent branches we can see the tree is splits by size where the right branch contains 47% of the data with a 92% probability of being Found. The final nodes use the geographical location to classify the remaining data where on the right branch containing 45% of the recordings the probability of being Found reaches 96%

Model Performance

To assess the performance of the model we will use the testing data: using the fitted decision tree, each observation in the testing data will be classified as either fall or found, and then we will check whether or not the classification was correct.

pred <-predict(mytree, test, type = 'class')
CM <- table(test$fall, pred)
CM
##        pred
##         Fell Found
##   Fell    87     4
##   Found   11    99
accuracy = sum(CM[1,1]+CM[2,2])/sum(CM)
precision = CM[2,2]/sum(CM[,2])
recall = CM[2,2]/sum(CM[2,])
specificity = CM[1,1]/sum(CM[1,])

accuracy
## [1] 0.9253731
precision
## [1] 0.961165
recall
## [1] 0.9
specificity
## [1] 0.956044

Conclusion

By balancing the data (the 50/50 split) we intentionally made the model less at predicting the reality but much better at telling you the why. Small meteorites (low mass) are much more likely to be destroyed by the environment (wind/water) or go unnoticed. Large meteorites are more resistant as they survive the environment longer, making them much more likely to be “Found” years or decades later. Also before 1973 there was a clear change in the way we look at meteorites, we stopped waiting for the meteorites just to fall and we decided to look for them. As we started to find more meteorites we also found new types that make our classification much wider than before this is why after 1973 the tree develops more leaves as this classification started to get more complex. A more meaningful way to interpret this tree is that the model has detected the exact moment when we stopped being just mere observers and decide to actively search for what falls from the space