# Cubic Clustering Criterion in Python

Does anybody know if any package calculates the Cubic Clustering Criterion (CCC) index and the Approximate Expected R Square (http://documentation.sas.com/?docsetId=emref&docsetTarget=n1dm4owbc3ka5jn11yjkod7ov1va.htm&docsetVersion=14.3&locale=en) in Python or have any example of the manual calculation of these metrics?

• Unfortunalely I don't know, for Python. But I once programmed CCC for SPSS (just FYI). – ttnphns Feb 8 at 13:26
• Thank you @ttnphns! Just found the solution. – Jony Zambrano Feb 9 at 8:12

I managed to create a function to calculate the CCC and Expected R Squared using this version in R - Cubic Clustering Criterion using R update. The input is a Dataframe with the features and the cluster assigned to each observation (as integer):

def ccc_calculation(c_r:"Dataframe with features and cluster number in the last column"):
n_c = int(c_r.iloc[:,-1:].max() + 1) # Number of clusters - assuming that there is a zero cluster - qq
n_o = c_r.shape[0] # Number of observations - nn
n_f = c_r.shape[1] - 1 # Number of clusters - pp
m_f = c_r.iloc[:,:-1].to_numpy() # Features matrix - jeu
m_c = c_r.iloc[:,-1:].to_numpy() # Clusters matrix
t_t = np.matmul(m_f.transpose(), m_f) # tt - p
v_e, m_e = np.linalg.eig(t_t/(n_o-1)) # Eigen values and eigen matrix
s_e = np.sqrt(v_e) # Square of eigen values
s_s = s_e[s_e != 0] # Exclude zeros - ss
v_v = np.prod(s_s) # Product of eigen values
m_z = sk.OneHotEncoder(sparse=False).fit_transform(m_c) # Z matrix with hot encoded clusters
m_x = np.matmul(np.matmul(np.linalg.inv(np.matmul(m_z.transpose(), m_z)), m_z.transpose()), m_f) # Matrix x
m_b = np.matmul(np.matmul(np.matmul(m_x.transpose(), m_z.transpose()), m_z), m_x) # Matrix b
m_w = np.subtract(t_t,m_b) # Matrix w
r_2 = 1 - (np.sum(np.diag(m_w))/np.sum(np.diag(t_t))) # R squared
m_u = s_s/(v_v/(n_c))**(1/n_f) # Matrix u
v_p = min(len(m_u[m_u>=1]), n_c - 1) # Value p
if v_p > 0 & v_p < n_f:
v_a = np.prod(s_s[0:v_p]) # Product of array
m_u = s_s/(v_a/(n_c))**(1/v_p) # Adjusted matrix u
b_1 = np.sum(1/(m_u[0:v_p] + n_o)) # b1
b_2 = np.sum((m_u[v_p:n_f]**2)/(m_u[v_p:n_f] + n_o)) # b2
e_r_2 = 1-((b_1+b_2)/np.sum(m_u**2))*(((n_o-n_c)**2)/n_o)*(1+4/n_o) # Expected r squared
c_c_c = mt.log((1-e_r_2)/(1-r_2))*(mt.sqrt(n_o*v_p/2)/((0.001+e_r_2)**1.2)) # CCC
else:
b_1 = np.sum(1/(m_u + n_o)) # b1
e_r_2 = 1-(b_1/sum(m_u**2))*(((n_o-n_c)**2)/n_o)*(1+4/n_o) # Expected r squared
c_c_c = mt.log((1-e_r_2)/(1-r_2))*(mt.sqrt(n_o*n_f/2)/((0.001+e_r_2)**1.2)) # CCC
return e_r_2, c_c_c