Correlations in the Measles Dataset vis-à-vis Practice assignment.

[Excitedwanderer]

Hello everyone, the intent behind starting this Q&A is to hopefully receive some constructive criticism or input regarding my (R code) solution to the measles practice assignment, and to also give a solution (and verify the correctness of it) to the bonus question therein.

I have commented my code and I will give a basic rundown of what its doing and how it works, my intent while writing the solution was to avoid loops and take advantage of inbuilt vectorized operations and functions of R.

The code:

# Store the dataframe in the variable measles
measles <- read.csv("https://deepayan.github.io/BSDS/2025-01-DARP/slides/data/measles.csv")

# Apply log transform to the values in the rate column
measles$rate<- log(measles$rate)

# Get unique state names and store it in the vector states.
states <- unique(measles$state)

# Get all unique pairs of states and store it as a matrix in UniqStatePairs.
UniqStatePairs <- t(combn(states,2))

# For any given pair of states, ComputeCorr finds the correlation between their rates.
ComputeCorr <- function(pair){
  state1 <- pair[1]
  state2 <- pair[2]
  state1df <- subset(measles, state == state1)
  state2df <- subset(measles, state == state2)
  
  mergedf <-merge(state1df,state2df, by="year", suffixes= c("_A","_B"))
  
  cor(mergedf$rate_A,mergedf$rate_B,use="complete.obs")
}

# Apply the function ComputeCorr to all unique pairs of states
Corr_rates <-apply(UniqStatePairs,1,ComputeCorr)


UniqStatePairs <- as.data.frame(UniqStatePairs)  # Convert the matrix to data frame
colnames(UniqStatePairs) <- c("State1", "State2")  # Rename the columns
UniqStatePairs$Correlation <- Corr_rates # Add correlation to each row

# Find the rows with minimum and maximum correlation
max_corr_row <- UniqStatePairs[which.max(UniqStatePairs$Correlation), ]
min_corr_row <- UniqStatePairs[which.min(UniqStatePairs$Correlation), ]

print(max_corr_row) 
print(min_corr_row) 

A couple of points where some explanation is required:
What does t(combn(states,2)) do? combn(states,2) generates all possible unordered pairs of states as a matrix where each column contains a pair of states, then the t function takes the transpose of this matrix so that the resulting matrix has in its each row a pair of states.

The ComputeCorr function takes as input a vector of length 2 with state names as its entries and returns the correlation between the rates of the two states. It does this by merging, by the year column, the subsets in the original dataset that contain entries of only those two states. Note further that use ="complete.obs" drops any rows containing NA values so that we get an appropriate correlation value (and not just NA).

We apply the function ComputeCorr rowwise to UniqStatePairs and assign the resulting vector of correlation values to Corr_rates. Afterwards we just convert UniqStatePairs to a dataframe and add Corr_rates as a new column. This gives us a dataframe with all unordered state pairs and their correlation values. After this getting the maximum and minimum correlation state pairs is straightforward.

Solution to the Bonus Question

Before thinking of the solution, I wanted to first understand what problem does the question pose. Supposing that we have removed all rows containing NA as a rate value, then what obstacle does this pose for us if we want to calculate the correlation between two states' rates? Clearly, what might happen is that one state may contain an NA value for the year 1996 while the other for 1992, so we don't have the paired data for those two years and we have to exclude those while calculating correlation. At least this is what I interpreted the problem to be, I would like to know if I may be missing something.

Assuming my interpretation is correct, we do not need to modify the given code at all, and in fact it returns the same output for the measles dataframe with NA values removed. To understand why this happens we take a closer look at ComputeCorr. Within the function when we merge the two subsets of the measles dataframe, we do it by = "year" and this ensures that the resulting merged dataframe has only those years in the year column which exists in both subset dataframes (basically an inner join). We get the solution this way to the bonus question without any extra effort! (contingent on my understanding being correct)

I would deeply appreciate any corrections, optimisations or simplifications in this approach, and it would be equally as wonderful to see alternate approaches as well. More than anything I would like to know whether my solution to the bonus question is correct.

[Excitedwanderer]

@deepayan Professor, does this satisfactorily answer the bonus question?

[deepayan]

Yes, this is fine. I would suggest using expand.grid() instead of combn() just to avoid any mistakes with the ordering of the combinations.

ustates <- unique(measles$state)
StatePairs <- expand.grid(State1 = ustates, State2 = ustates)
StatePairs$Correlation <- apply(StatePairs, 1, ComputeCorr)

An alternative approach is the following, which takes advantage of the fact that cor() works on a data.frame or matrix, giving the pairwise correlation matrix.

## option 1: Use the fact that all state-year combinations are
## represented in a regular order. This would not work when the NA
## rows are removed altogether from the dataset, or even if the rows 
## are scrambled.

measlesMatrix1 <- matrix(measles$rate, nrow = 76, ncol = 51, byrow = TRUE)
cor1 <- cor(measlesMatrix1, use = "pairwise.complete.obs")
dimnames(cor1) <- list(ustates, ustates)

## options 2: use xtabs(), which can handle missing combinations etc.

measlesMatrix2 <- xtabs(rate ~ year + state, measles)
cor2 <- cor(measlesMatrix2, use = "pairwise.complete.obs")

## check that these give the same results

all.equal(cor1, cor2)

## need to ignore diagonal entries when finding maximum

diag(cor1) <- NA
array2DF(cor1)[which.max(cor1), ]
array2DF(cor1)[which.min(cor1), ]

An interesting follow up problem would be to continue this process
by finding the next pair of 'most similar states', treating the
first pair found as a single unit. This would need some way to
compare one group of (possibly multiple) states with another group.

A very crude way of doing this is the following

d <- 1 - cor2 ## d is a 'dissimilarity matrix'
h <- hclust(as.dist(d), method = "complete")
plot(h)