Putting the Elements of a Cluster into a Dataframe

Below I use functional principal component analysis to reduce the dimensionality of a pre-1989 time series. After, I cluster over the fPCA scores. Consider the following, which sets up the problem:

library(splines)
library(fda)
library(dplyr)
set.seed(1420)
data <-read.csv("https://raw.githubusercontent.com/synth-inference/synthdid/master/data/california_prop99.csv",
header=TRUE, sep=";")
# Imports our dataset
data = subset(data, select = c(State, Year,PacksPerCapita,treated))
# Restricts our dataset to the outcome, time, and panel
preperiod = nrow(data[data$treated== 1 & data$State=="California", ])
# Now we count the number of T_0 periods...
data <-data %>% select(-(treated))
#The treatment variable is no longer needed.
nptperiod=preperiod+2
if(preperiod < 3){
# If we have less than 3 pre-periods, stop.
stop("error message to print")   
}
# Step 0.2: Rehshape our dataframe to wide.
trainframe = reshape(data, idvar = "State", timevar = "Year", direction = "wide")
# Step 1: Define X as the pre-intervention
#period for both treated unit and
#non-treated units in the data set
X = trainframe[,2:nptperiod]
last = nptperiod-1
x = seq(0,1,length=last)
# Step 1.01: Define our spline for the eigenfunction
splinebasis = create.bspline.basis(c(0,1),10)
smooth = smooth.basis(x,t(X),splinebasis)
# Step 1.2: Compute FPC scores xi for all units
# that explain most of the variation in the data
Xfun = smooth$fd
#Do PCA over it
pca = pca.fd(Xfun, 10)
var.pca = cumsum(pca$varprop)
nharm = sum(var.pca < 0.95) + 1
pc = pca.fd(Xfun, nharm)
# Step 2: Apply K-Means algorithm on
# FPC-scores and extract the donor pool
#Make our fPCA scores as matrices
cluster_x = as.matrix(pc$scores)
cluster_x <- scale(cluster_x)
# 2.1: Calculating the Silhouette statistics
library(cluster)
k.max = 8
sil = rep(0, k.max)
# Compute the average silhouette width 
for(i in 2:k.max){
tmp = kmeans(cluster_x, centers = i, nstart = 10)
ss <- silhouette(tmp$cluster, dist(cluster_x))
# Our Silhouette stat
sil[i] <- mean(ss[, 3])
}
# Now, we get the row of the vector
# that has the largest Silhouette stat
optnum = which.max(sil)
# 2.2: K-Means Cluster Analysis
fit <- kmeans(cluster_x, optnum) # 5 cluster solution
# 5 cluster solution?
trainframe <- data.frame(trainframe, fit$cluster)
cols=ncol(trainframe)
colm1 = cols-1

The cluster California (the treated unit) belongs in is cluster 1, and I want to keep states which only are in cluster 1. But, when I check to see if it's in new_data after clustering like so

trainframe[,c(1,cols)]
treat_cluster = fit$cluster[11]
new_data = trainframe[trainframe[,cols]==treat_cluster,1:colm1]
istr= sum(str_detect(new_data$fullname, '^California$')) > 0
if(istr == "FALSE"){
# If California isn't there, stop.
stop("You MUST have California in the cluster.")   
}

R tells me that California isn't found, which it should be. How do I save cluster 1 into its own dataframe? I suspect the issue has to do with treat_cluster = fit$cluster[11].

As promised, here's the solution. We begin with the setup

library(splines)
library(fda)
library(ggplot2)
library(dplyr)
library(stringr)
set.seed(1420)
data <-read.csv("https://raw.githubusercontent.com/synth-inference/synthdid/master/data/california_prop99.csv",
header=TRUE, sep=";")
# Imports our dataset
data = subset(data, select = c(State, Year,PacksPerCapita,treated))
# Restricts our dataset to the outcome, time, and panel
preperiod = nrow(data[data$treated== 0 & data$State=="California", ])
# Now we count the number of T_0 periods...
data <-data %>% select(-(treated))
#The treatment variable is no longer needed.
nptperiod=preperiod+2
if(preperiod < 3){
# If we have less than 3 pre-periods, stop.
stop("error message to print")   
}
# Step 0.2: Rehshape our dataframe to wide.
trainframe = reshape(data, idvar = "State", timevar = "Year", direction = "wide")
# Step 1: Define X as the pre-intervention
#period for both treated unit and
#non-treated units in the data set
X = trainframe[,2:nptperiod]
last = nptperiod-1
x = seq(0,1,length=last)
# Step 1.01: Define our spline for the eigenfunction
splinebasis = create.bspline.basis(c(0,1),10)
smooth = smooth.basis(x,t(X),splinebasis)
# Step 1.2: Compute FPC scores xi for all units
# that explain most of the variation in the data
Xfun = smooth$fd
#Do PCA over it
pca = pca.fd(Xfun, 10)
var.pca = cumsum(pca$varprop)
nharm = sum(var.pca < 0.95) + 1
pc = pca.fd(Xfun, nharm)

where we estimate the functional principal components. Next we k-means cluster over them.

# Step 2: Apply K-Means algorithm on
# FPC-scores and extract the donor pool
#Make our fPCA scores as matrices
cluster_x = as.matrix(pc$scores)
cluster_x <- scale(cluster_x)
# 2.1: Calculating the Silhouette statistics
library(cluster)
k.max = 8
sil = rep(0, k.max)
# Compute the average silhouette width 
for(i in 2:k.max){
tmp = kmeans(cluster_x, centers = i, nstart = 10)
ss <- silhouette(tmp$cluster, dist(cluster_x))
# Our Silhouette stat
sil[i] <- mean(ss[, 3])
}
# Now, we get the row of the vector
# that has the largest Silhouette stat
optnum = which.max(sil)
# 2.2: K-Means Cluster Analysis
fit <- kmeans(cluster_x, optnum) # 5 cluster solution
trainframe <- data.frame(trainframe, fit$cluster)
cols=ncol(trainframe)
colm1 = cols-1
cluster_num = sum(trainframe[which(trainframe[[1]] =='California'), cols])

The last step here is key- I extract the row-value of the cluster column if column 1, our unit column, is California. Next, we match the cluster California falls within with its place in fit$cluster, the list of clusters in order.

trainframe[,c(1,cols)]
position = match(c(cluster_num),fit$cluster)
treat_cluster = fit$cluster[position]
new_data = trainframe[trainframe[,cols]==treat_cluster,1:colm1]

We now have a dataframe consisting of the optimal subset of donors.