#Article 2
# Factor Analysis
library(psych)
tetrachoric(Article2[c("Quality","Organization","Product","ProductImprovement","Efficiency","Flexibility")])
corr<-tetrachoric(Article2[c("Quality","Organization","Product","ProductImprovement","Efficiency","Flexibility")])$rho
cortest.bartlett(corr,1139)
KMO(corr)
fa.parallel(Article2[c("Quality","Organization","Product","ProductImprovement","Efficiency","Flexibility")],cor="tet")
fa <- fa(r = corr,nfactors = 2,n.obs = 1139, rotate = "varimax")
fa[["loadings"]]
fa[["communalities"]]
fs <- factor.scores(Article2[c("Quality","Organization","Product","ProductImprovement","Efficiency","Flexibility")], fa)
fs <- fs$scores
df <- cbind(Article2,fs)

# Alluvial diagram
library(ggplot2)
library(alluvial)
df$Freq <- 1
Aggregates1<- aggregate(Freq~ Quality + Organization + ProductImprovement + Product, data=df, FUN=sum)
Aggregates2<- aggregate(Freq~ Efficiency + Flexibility, data=df, FUN=sum)
par(mfrow=c(1,2))
alluvial(Aggregates1[, -ncol(Aggregates1)], freq = Aggregates1$Freq)
alluvial(Aggregates2[, -ncol(Aggregates2)], freq = Aggregates2$Freq, cw=0.03)

#Jitter plot
par(mfrow=c(1,2))
plot(jitter(df$`prod_num`),jitter(df$`factor1`),col="#000000",pch=16,xlab="Number of product-oriented innovations", ylab = "Factor 1")
plot(jitter(df$`proc_num`),jitter(df$`factor2`),col="#000000",pch=16,xlab="Number of production-oriented innovations", ylab = "Factor 2")

#Article 3
#China map
library(plyr)
library(ggplot2)
china_data <- join(china_map,mydata,type="full")
ggplot(china_data,aes(x=long,y=lat))+geom_polygon(aes(group=group,fill=Number),colour="grey",size=0.01)+scale_fill_gradient(low="yellow",high="red")+coord_map("polyconic")+geom_text(aes(x=jd,y=wd,label=NAME),data=province,colour="black",size=2.5)

#Cluster analysis on innovation activities
library(cluster)
library(NbClust)
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")
ac <- function(x) {agnes(Article3[c("af1","af2")], method = x)$ac}
sapply(m, ac)
NbClust(Article3[c("af1","af2")], distance = "euclidean", method = "ward.D2", max.nc = 10, index = "all")
d <- dist(Article3[c("af1","af2")], method = "euclidean")
final_clust <- hclust(d, method = "ward.D2" )
groups <- cutree(final_clust, k=3)
Article3 <- cbind(Article3, acluster = groups)

#Cluster analysis on innovation types
library(cluster)
library(NbClust)
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")
ac <- function(x) {agnes(Article3[c("tf1","tf2")], method = x)$ac}
sapply(m, ac)
NbClust(Article3[c("tf1","tf2")], distance = "euclidean", method = "ward.D2", max.nc = 10, index = "all")
d <- dist(Article3[c("tf1","tf2")], method = "euclidean")
final_clust <- hclust(d, method = "ward.D2" )
groups <- cutree(final_clust, k=3)
Article3 <- cbind(Article3, tcluster = groups)