library(ca) library(psych) data <- read.table("data/sciencedoctorates.txt",header=T,sep="\t",row.names=1) # View(data) dim(data) # Again, remove the total col/row SD <- data[-dim(data)[1],-dim(data)[2]] dim(SD) # To interpret correspondence analysis, the first step is to evaluate whether there is a significant # dependency between the rows and columns. n <- sum(SD) v1 <- matrix(colSums(SD),nrow=1) v2 <- matrix(rowSums(SD),ncol=1) #theoretical frequencies under independence E <- v2 %*% v1/n I <- dim(SD)[1] J <- dim(SD)[2] s <- 0 #chi-square statistic for(i in 1:I){ for(j in 1:J){ s <- s + ( SD[i,j]-E[i,j] )^2/(E[i,j]) } } pchisq(s,df=((I-1)*(J-1)),lower.tail=F) # H0: Discipline and Year are independent # H1: Discipline and Year are not independent chisq.test(SD) # there is evicende that there is statistically # significant association between the number # of doctors graduated and the year (in USA) #ca-function help(ca) SD.ca = ca(SD) #set nd=8 to make all columns visible in summary names(SD.ca) SD.ca\$colnames SD.ca\$rownames SD.ca\$N SD.ca\$sv #The square roots of singular values related to the PCA transformation for rows/cols # (how much variation explained by the principal components) # for symmetric matrices, singular values = |eigenvalues| # (here sv's are used since the package uses svd instead of eigen) ######################## # Inertia ######################## # How much of the total "variation" the specific variable explains # i.e. how much it contributes to the chi-squared statistic # Inertia is the chi squared statistic divided by n SD.ca\$rowinertia SD.ca\$colinertia # You can get the single row-inertia values by fixing i (the row index) s2 <- 0 #For physics, i =3 i <- 3 for(j in 1:J){ s2 <- s2 + ( SD[i,j]-E[i,j] )^2/(E[i,j]) } s2/n #Note that sum(SD.ca\$rowinertia) sum(SD.ca\$colinertia) #is the same as sum(SD.ca\$sv^2) s/n SD.ca\$rowinertia/sum(SD.ca\$sv^2) #these proportional values are the ones seen in summary(SD.ca) SD.ca\$colinertia/(s/n) ############################################### #ca manually ################################### #theoretical relative frequencies under independence Ef <- E/n #observed relative frequencies SDf = SD/n #the matrix Z Z = (SDf - Ef)/sqrt(Ef) class(Z) Z = as.matrix(Z) #matrices V and W V = t(Z)%*%Z W = Z%*%t(Z) #chi-squared test statistic n*tr(V) variances = eigen(V)\$values SD.ca\$sv^2 components = eigen(V)\$vectors #V and W has the same non-zero eigenvalues eigen(W)\$values components2 = eigen(W)\$vectors #forming the matrix R f1 = v1/n #SD.ca\$colmass f2 = v2/n #SD.ca\$rowmass one = matrix(rep(1,nrow(SD)), ncol=1) sifting = one%*%sqrt(f1) scaling = f2 %*% sqrt(f1) R = SDf/scaling - sifting #rowcoordinates rowcoord = as.matrix(R)%*%components #omit the dimension corresponding to zero eigenvalue rowcoord = rowcoord[,-8] summary(SD.ca) #standardized rowcoordinates stand = matrix(rep(1,nrow(SD), ncol=1))%*%matrix(variances[-8], nrow=1) standrowcoord = rowcoord/sqrt(stand) SD.ca\$rowcoord #forming the matrix C one2 = matrix(rep(1,ncol(SD)), nrow=1) sifting2 = sqrt(f2)%*%one2 scaling2 = sqrt(f2) %*% f1 C = SDf/scaling2 - sifting2 #colcoordinates colcoord = t(C)%*%components2 #omit dimensions corresponding to zero eigenvalues colcoord = colcoord[,-c(8:ncol(colcoord))] summary(SD.ca) #standardized colcoordinates stand2 = matrix(rep(1,ncol(SD), ncol=1))%*%matrix(variances[-8], nrow=1) standcolcoord = colcoord/sqrt(stand2) SD.ca\$colcoord ## The chi-squared distances from the "center", where variables close to center do not deviate from the # independence assumption. sqrt(rowSums(rowcoord^2)) SD.ca\$rowdist #coldistances sqrt(rowSums(colcoord^2)) SD.ca\$coldist #################################################### SD.ca summary(SD.ca) # Note that the quantities are multiplied by 1000 # Quality of representation = as in the lecture slides, but here we consider the angle between profiles and # the plane spanned by the two first principal components # Squared correlations = quality of representation from lecture slides # also, the sum of the squared correlations is the quality of representation. # ctr = contribution in forming that ca-component (contributions sum to 1) # important variables related to forming the specific component have a high ctr # k=1 and k=2 are the coordinates on the plot names(summary(SD.ca)\$rows) # Contribution of engineering to the second axis f2[1]*SD.ca\$rowcoord[1,2]^2 #recall that these coordinates are already scaled # If the rows and columns were independent, ctr would be same for every variable # Squared correlation of biology with the second component d2 = rowcoord[6,]^2 d2[2]/sum(d2) #Note that the following plots unscaled coordinates (principal coordinates) and hence, deduction based on the plot is questionable plot(SD.ca,arrows=c(T,T),map="symmetric") #plot(SD.ca,arrows=c(T,T),map="symmetric",dim=c(2,4)) #Instead try e.g. following commands plot(SD.ca, arrows=c(T,T), map="rowprincipal") #standard column coordinates plot(SD.ca, arrows=c(T,T), map="rowgreen") #standard column coordinates scaled with square root of the column masses plot(SD.ca, arrows=c(T,T), map="colgab", dim=c(3,4)) #standard row coordinates scaled with the row masses plot(SD.ca, arrows=c(T,T), map="rowgab") # If two row-variables are close on the picture, they have a similar profile, # the same is true for column-variables # Distant row/column-variables have different profiles # Variables distant from the origin represent variables different from the average profile # these are usually the most interesting ones # Now you can again try to interpret the dimensions. # 1st dim splits the sciences into soft/hard # 2nd dim splits the sciences into more formula heavy(math,physics,engineering) vs # the more experimental ones (chemistry,agriculture,earthsciences) # same for different years # You can also try this: plot3d.ca(SD.ca)