library(likert)
library(plot.matrix)
library(semPlot)
library(semTools)
library(MASS)

library(fields)

# x is a matrix containing the data
# method : correlation method. "pearson"" or "spearman"" is supported
# removeTriangle : remove upper or lower triangle
# results :  if "html" or "latex"
# the results will be displayed in html or latex format
corstars <-function(x, method=c("pearson", "spearman"), removeTriangle=c("upper", "lower"),
                    result=c("none", "html", "latex")){
  #Compute correlation matrix
  require(Hmisc)
  x <- as.matrix(x)
  correlation_matrix<-rcorr(x, type=method[1])
  R <- correlation_matrix$r # Matrix of correlation coeficients
  p <- correlation_matrix$P # Matrix of p-value 
  
  ## Define notions for significance levels; spacing is important.
  #mystars <- ifelse(p < .0001, "****", ifelse(p < .001, "*** ", ifelse(p < .01, "**  ", ifelse(p < .05, "*   ", "    "))))
  mystars <- ifelse(p < .01, "**  ", ifelse(p < .05, "*   ", "    "))
  
  ## trunctuate the correlation matrix to two decimal
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1]
  
  ## build a new matrix that includes the correlations with their apropriate stars
  Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x))
  diag(Rnew) <- paste(diag(R), " ", sep="")
  rownames(Rnew) <- colnames(x)
  colnames(Rnew) <- paste(colnames(x), "", sep="")
  
  ## remove upper triangle of correlation matrix
  if(removeTriangle[1]=="upper"){
    Rnew <- as.matrix(Rnew)
    Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
    Rnew <- as.data.frame(Rnew)
  }
  
  ## remove lower triangle of correlation matrix
  else if(removeTriangle[1]=="lower"){
    Rnew <- as.matrix(Rnew)
    Rnew[lower.tri(Rnew, diag = TRUE)] <- ""
    Rnew <- as.data.frame(Rnew)
  }
  
  ## remove last column and return the correlation matrix
  Rnew <- cbind(Rnew[1:length(Rnew)-1])
  if (result[1]=="none") return(Rnew)
  else{
    if(result[1]=="html") print(xtable(Rnew), type="html")
    else print(xtable(Rnew), type="latex") 
  }
} 

oldnames = paste0("P",1:8)
newnames<- structure(c("PCK","TPK","TK","TCK","TPACK","CK","PK","BNCL")
                     ,names=oldnames )

dat <- read.csv("Primary school language and literacy teachers' knowledge in teaching literacies associated with a ra_August 10, 2020_15.38 - Sheet1.csv",skip=1,as.is=T)

labs <- c("Strongly disagree","Disagree","Somewhat disagree","Neither agree nor disagree"
          ,"Somewhat agree","Agree","Strongly agree") 

names(dat) <- gsub("\\."," ",names(dat))
names(dat) <- gsub(" Mark only one circle ","",names(dat))

for (i in 1:dim(dat)[2]) {
  col = dat[,i]
  if(sum(col %in% labs)>1) {
     dat[,i] <- factor(col,levels=labs,ordered=T)
  }
}


dat$age <- factor(dat$`Age  `)


ocols <- c()
for (i in 1:8) {
  part = paste0("P",i)
  print(part)
  cols <- grep(part,names(dat))
  ocols <- c(ocols,cols)
  png(paste0(part,".png"),width=800,height=400)
  print(plot(likert(dat[,cols]),type="density",title="MAIN"))
  dev.off()
}
  
png("Fig_1-8.png",height=400,width=800)
layout(mat=matrix(1:8,nrow=4))
ocols <- c()
for (i in 1:8) {
  part = paste0("P",i)
  print(part)
  cols <- grep(part,names(dat))
  ocols <- c(ocols,cols)
  print(plot(likert(dat[,cols]),type="density",title="MAIN"))
}

dev.off()



dat$exp <- dat$`How long have you been teaching in a primary school in your teaching career  Provide your answers in this format   XX years XX months  ` 

dat$exp <- toupper(dat$exp)
dat$expy <- sapply(strsplit(dat$exp,"Y"),function(l) l[1])
dat$expm <- sapply(strsplit(dat$exp,"Y"),function(l) l[2])

dat$expy <- gsub("~","",dat$expy)
dat$expy <- gsub("\\+","",dat$expy)
dat$expy <- ifelse(dat$expy=="2806","28.5",dat$expy)
dat$expy <- ifelse(dat$expy=="0306","3.5",dat$expy)
dat$expy <- ifelse(grepl("M",dat$expy),
                   as.numeric(gsub("MONTHS","",dat$expy))/12,
                   dat$expy )

dat$expm <- gsub("[^0-9.]","",dat$expm)
dat$expm <- ifelse(is.na(dat$expm)|dat$expm=="",0,dat$expm)
dat$expm <- as.numeric(dat$expm)/12
dat$expy <- as.numeric(dat$expy)+dat$expm

dat[,c(63,64,65)]

dat$permteach <- dat[,4] == "In-service teacher (permanent)"
dat$casual <- dat[,4] =="Casual teacher (contracted for no more than a teaching term)"

dat$NSW <- FALSE
dat$NSW <- ifelse(grepl("New",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("NSW",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("NDW",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("Nsw",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("nsw",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("Marrickville",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("Granville",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(grepl("Tallowood",dat[,6]),TRUE,dat$NSW)
dat$NSW <- ifelse(dat[,6]=="N/A" | dat[,6]=="",NA,dat$NSW)



### CFA ###
dataq <- dat[,9:53]
names(dataq) <- substr(names(dataq),1,4)
names(dataq) <- gsub(" ","q",names(dataq))
library(lavaan)
library(lavaanPlot)

dataq <- sapply(dataq,as.numeric)

m1 <- 'PCK =~ P1q1 + P1q2 + P1q3 + P1q4 + P1q5 + P1q6
        TPK =~ P2q1 + P2q2 + P2q3 + P2q4 + P2q5 + P2q6
        TK =~ P3q1 + P3q2 + P3q3 + P3q4
        TCK =~ P4q1 + P4q2 + P4q3 + P4q4
        TPACK =~ P5q1 + P5q2 + P5q3 + P5q4 + P5q5 + P5q6
        CK =~ P6q1 + P6q2 + P6q3 + P6q4 + P6q5 + P6q6
        PK =~ P7q1 + P7q2 + P7q3 + P7q4 + P7q5 + P7q6
        BNCL =~ P8q1 + P8q2 + P8q3 + P8q4 + P8q5 + P8q6 + P8q7
        '
model1 <- cfa(m1,data=dataq)
summary(model1,fit.measures=T)

#r=50
#png("cfa.png",width=480*r,height=960*r)
#semPaths(model1,what="est", title = FALSE,layout="circle")
#dev.off()

write.csv(round(t(reliability(model1))[,c(1,5)],2),"reliability1.csv")

rownames(dataq) <- 1:nrow(dataq)
dataqsub <- dataq[rowSums(is.na(dataq))==0,]
efa=factanal(dataqsub,factors=8,rotation="varimax")
efa

# Drop TCK q4

m2 <- 'PCK =~ P1q1 + P1q2 + P1q3 + P1q4 + P1q5 + P1q6
        TPK =~ P2q1 + P2q2 + P2q3 + P2q4 + P2q5 + P2q6
        TK =~ P3q1 + P3q2 + P3q3 + P3q4
        TCK =~ P4q1 + P4q2 + P4q3 
        TPACK =~ P5q1 + P5q2 + P5q3 + P5q4 + P5q5 + P5q6
        CK =~ P6q1 + P6q2 + P6q3 + P6q4 + P6q5 + P6q6
        PK =~ P7q1 + P7q2 + P7q3 + P7q4 + P7q5 + P7q6
        BNCL =~ P8q1 + P8q2 + P8q3 + P8q4 + P8q5 + P8q6 + P8q7
        '
model2 <- cfa(m2,data=dataq)
summary(model2,fit.measures=T)

write.csv(round(t(reliability(model2))[,c(1,5)],2),"reliability2.csv")


png("Correlation.png")
r <- cor(dataq,method="kendall",use="pairwise.complete.obs")
col2 <- colorRampPalette(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                          "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                          "#4393C3", "#2166AC", "#053061"))

image.plot(r,axes=F,col=col2(200),breaks=seq(-1,1,length.out=201))
box()
nq = c(0,6,6,4,4,6,6,6,7)
b=sum(nq)
tl=1+1/(b-1)
ln = tl*cumsum(nq/sum(nq))-0.5/(b-1)
abline(v=ln)
abline(h=ln)
mtext("Correlation between items",cex=1.2,line=1)
axis(1,at=ln[-9]+diff(ln)/2,labels=newnames,las=2)
axis(2,at=ln[-9]+diff(ln)/2,labels=newnames,las=2)
dev.off()

library(ggcorrplot)
png("Correlation_factors.png")
r <- cor(newdat[,oldnames])
row.names(r) <- newnames
colnames(r) <- newnames
ggcorrplot(r)
dev.off()

library(corrplot)

oldcolnames <- colnames(dataq)
for(i in oldnames) {
  colnames(dataq) <- gsub(i,newnames[i],colnames(dataq))
}

rstar <- corstars(newdat[,oldnames])
for (i in 1:length(rstar)) rstar[i,i]="1"
row.names(rstar) = newnames[oldnames]
names(rstar) = newnames[oldnames][-8]

write.csv(rstar,"cortable.csv")

rall <- cor(dataq,method="kendall",use="pairwise.complete.obs")
write.csv(round(rall,2),"allpairwise_corrs.csv")

#CK and PK
sub=c("CK","PK")
oldsub = names(newnames)[newnames %in% sub]
cols <- grepl(oldsub[1],oldcolnames) | grepl(oldsub[2],oldcolnames)

r1 <- cor(dataq[,cols],method="kendall",use="pairwise.complete.obs")
png("corrplot_CK_PK.png",height=720,width=720)
corrplot.mixed(r1,tl.srt=45,tl.cex=0.8)
segments(0.5,6.5,12.5,6.5,xpd=T)
segments(6.5,0.5,6.5,12.5,xpd=T)
dev.off()


#TPK and TPACK
sub=c("TPK","TPACK")
oldsub = names(newnames)[newnames %in% sub]
cols <- grepl(oldsub[1],oldcolnames) | grepl(oldsub[2],oldcolnames)

r2 <- cor(dataq[,cols],method="kendall",use="pairwise.complete.obs")
png("corrplot_TPK_TPACK.png",height=720,width=720)
corrplot.mixed(r2,tl.srt=45,tl.cex=0.8)
segments(0.5,6.5,12.5,6.5,xpd=T)
segments(6.5,0.5,6.5,12.5,xpd=T)
dev.off()


# CK, PCK and TCK
sub=c("CK","TCK","PCK")
oldsub = names(newnames)[newnames %in% sub]
cols <- grepl(oldsub[1],oldcolnames) | grepl(oldsub[2],oldcolnames)| grepl(oldsub[3],oldcolnames)

r3 <- cor(dataq[,cols],method="kendall",use="pairwise.complete.obs")
png("corrplot_CK_TCK_PCK.png",height=720,width=720)
corrplot.mixed(r3,tl.srt=45,tl.cex=0.8)
segments(0.5,c(6.5,10.5),16.5,c(6.5,10.5),xpd=T)
segments(c(6.5,10.5),0.5,c(6.5,10.5),16.5,xpd=T)
dev.off()




# Drop q4
dim(dat)
dat <- dat[,!grepl("P4 4",names(dat))]
dim(dat)

dim(dataq)
dataq <- dataq[,colnames(dataq) != "P4q4"]
dim(dataq)

newdat <- dat[,c("Gender  ","expy","Are you a current PETAA member ","age","casual","NSW")]
for (i in 1:8) {
  item = paste0("P",i)
  newdat[,item] <- rowMeans(as.data.frame(lapply(dat[,grep(item,names(dat))],as.numeric)),na.rm=T)
}




res <- data.frame(names=oldnames,
                  means=unlist(lapply(newdat[,oldnames],function(c) mean(c,na.rm=T))))
res$CI_L <- sapply(oldnames
                   ,function(c) confint(lm(data=newdat,as.formula(paste0(c,"~1"))))[1] )
res$CI_U <- sapply(oldnames
                   ,function(c) confint(lm(data=newdat,as.formula(paste0(c,"~1"))))[2] )
res <- res[order(res$means),]





pdf("CompareDomains.pdf")
par(mar=c(12,7,1,1))
plot(res$means,1:8,xlim=c(1,7),axes=F
     ,xlab="",ylab="",pch=19)
abline(v=1:7,col="gray")
abline(h=1:8,col="gray")

mtext("Domain",2,5,cex=1.5)
box()
axis(1,1:7,labs,las=2)
axis(2,1:8,paste0(res$names,":",newnames[res$names]),las=2)
arrows(res$CI_L,1:8,res$CI_U,1:8,code=3,angle=90,length=0.04,lwd=2)
dev.off()

newdat$age_n <- as.numeric(newdat$age)

names(newdat)[grep("PETAA",names(newdat))] <- "PETAA"  
newdat$Gender <- factor(newdat$`Gender  `)
newdat$Gender[newdat$Gender==""] <- NA
mfit <- manova(cbind(P1,P2,P3,P4,P5,P6,P7,P8)~age_n+expy+Gender+PETAA+casual+NSW,data=newdat)
summary(mfit)
summary.aov(mfit)
df <- summary(mfit)$stats
write.csv(round(df,3),"manova_results.csv")




layout(mat=matrix(1:16,nrow=8,byrow=T),heights=c(rep(1,7),1.5))
par(mar=c(0.1,4,0.1,0.1))
for(i in 1:8) {
  if (i==8) {
    par(mar=c(6,4,0.1,0.1))
  }
  f = as.formula(paste0(oldnames[i],"~expy+age_n"))
  fiti = lm(f,data=newdat)
  plot(newdat$expy,newdat[,oldnames[i]],xlab="Years of experience",ylab=
         paste0(oldnames[i],":",newnames[oldnames[i]]),ylim=c(1,7),xaxt="n")
  if(i==8) {
    axis(1)
  }
  pv = summary(fiti)$coefficients[2,4]
  b= summary(fiti)$coefficients[2,1]
  abline(fiti$coef[1]+mean(newdat$age_n)*fiti$coef[3],fiti$coef[2]
         ,col=ifelse(pv<(0.05),"black","gray"))
  legend("bottomright",legend=
    paste0("slope=",round(b,3)," p-value=",round(pv,3)))
  
  plot(newdat$age_n,newdat[,oldnames[i]],xlab="Age",ylab=
         paste0(oldnames[i],":",newnames[oldnames[i]]),ylim=c(1,7),xaxt="n")
  if(i==8) {
    axis(1,at=1:nlevels(dat$age),levels(dat$age),las=2)
  }
  pv = summary(fiti)$coefficients[3,4]
  b= summary(fiti)$coefficients[3,1]
  abline(fiti$coef[1]+mean(newdat$expy,na.rm=T)*fiti$coef[2],fiti$coef[3]
         ,col=ifelse(pv<(0.05),"black","gray"))
  legend("bottomright",legend=
           paste0("slope=",round(b,3)," p-value=",round(pv,3)))
}

dev.off()


newitems <- colnames(dataq)
for (n in oldnames) {
  newitems <- gsub(n,newnames[n],newitems)
}
colnames(dataqsub) <- newitems

expfit=lm(expy+0.5~P7+age_n,data=newdat)
imp <- function(row) {
  if(is.na(newdat$expy[row])) {
    return( max(0,predict(expfit,newdata=newdat[row,])) )
  } else {
    return( newdat$expy[row] )
  }
} 
newdat$expy_i <- sapply(1:nrow(newdat),imp)


forcl <- newdat[,c(oldnames,"age_n","expy_i")]

library(clValid)

test=clValid(forcl,2:4,clMethods=c("hierarchical","kmeans","pam"),validation="internal")
summary(test)
plot(test)


png("Internal Valdiation.png")
op <- par(no.readonly=TRUE)
par(mfrow=c(2,2),mar=c(4,4,3,1)) 
plot(test, legend=FALSE)
par(op)
dev.off()

kc <- kmeans(newdat[,c(oldnames)],centers=2)
kc

agekey <- data.frame(agel=levels(newdat$age),age=c(21,27.5,33,38,43,48,53,58,63,67))
newdat$agey <- agekey$age[match(newdat$age,agekey$agel)]
vars=c(oldnames,"agey","expy")
cl1 <- newdat[kc$cluster==1,]
res1 <- data.frame(names=vars,
                  means=unlist(lapply(cl1[,vars],function(c) mean(c,na.rm=T))))
res1$CI_L <- sapply(vars
                   ,function(c) confint(lm(data=cl1,as.formula(paste0(c,"~1"))))[1] )
res1$CI_U <- sapply(vars
                   ,function(c) confint(lm(data=cl1,as.formula(paste0(c,"~1"))))[2] )


res1_ae <- res1[9:10,]
res1 <- res1[1:8,]

cl2 <- newdat[kc$cluster==2,]
res2 <- data.frame(names=vars,
                   means=unlist(lapply(cl2[,vars],function(c) mean(c,na.rm=T))))
res2$CI_L <- sapply(vars
                    ,function(c) confint(lm(data=cl2,as.formula(paste0(c,"~1"))))[1] )
res2$CI_U <- sapply(vars
                    ,function(c) confint(lm(data=cl2,as.formula(paste0(c,"~1"))))[2] )
res2_ae <- res2[9:10,]
res2 <- res2[1:8,]

diff = res1$means-res2$means
res1 <- res1[order(diff),]
res2 <- res2[order(diff),]



pdf("CompareClusters.pdf")

par(mar=c(12,7,1,1))
plot(res1$means,c(1:8)*2,xlim=c(1,7),ylim=c(1,16),axes=F
     ,xlab="",ylab="",pch=19,col=c("black"))
points(res2$means,c(1:8)*2-1,pch=19,col=c("gray"))

abline(v=1:7,col="light gray")
abline(h=2*c(1:8)+0.5,col="light gray")

mtext("Variable",2,5,cex=1.5)
box()
axis(1,1:7,labs,las=2)
axis(2,c(1:8)*2-.5,paste0(res1$names,":",newnames[res1$names]),las=2)
arrows(res1$CI_L,2*c(1:8),res1$CI_U,2*c(1:8),code=3,angle=90,length=0.04,lwd=2)
arrows(res2$CI_L,2*c(1:8)-1,res2$CI_U,2*c(1:8)-1,code=3,angle=90,length=0.04,lwd=2,col="grey")
legend("bottomleft",col=c("black","grey"),pch=19,legend=c("cluster 1","cluster 2"))
dev.off()


pdf("CompareClusters2.pdf")

par(mar=c(12,7,1,1))
plot(res1_ae$means,c(1:2)*2,ylim=c(1,4),yaxt="n"
     ,xlim=range(res1_ae[,-1],res2_ae[,-1])
     ,xlab="",ylab="",pch=19,col=c("black"))
points(res2_ae$means,c(1:2)*2-1,pch=19,col=c("gray"))

#abline(v=1:7,col="light gray")
abline(h=2*c(1:2)+0.5,col="light gray")

mtext("Variable",2,5,cex=1.5)
box()
#axis(1,1:7,labs,las=2)
axis(2,c(1:2)*2-.5,res1_ae$names,las=2)
arrows(res1_ae$CI_L,2*c(1:2),res1_ae$CI_U,2*c(1:2),code=3,angle=90,length=0.04,lwd=2)
arrows(res2_ae$CI_L,2*c(1:2)-1,res2_ae$CI_U,2*c(1:2)-1,code=3,angle=90,length=0.04,lwd=2,col="grey")
legend("bottomleft",col=c("black","grey"),pch=19,legend=c("cluster 1","cluster 2"))
dev.off()


# AS TABLE!
vars=c("agey","expy",oldnames)

cl1 <- newdat[kc$cluster==1,]
res1 <- data.frame(M=unlist(lapply(cl1[,vars],function(c) mean(c,na.rm=T))))
res1$CI_L <- sapply(vars
                    ,function(c) confint(lm(data=cl1,as.formula(paste0(c,"~1"))))[1] )
res1$CI_U <- sapply(vars
                    ,function(c) confint(lm(data=cl1,as.formula(paste0(c,"~1"))))[2] )
res1$SD = unlist(lapply(cl1[,vars],function(c) sd(c,na.rm=T)))
row.names(res1)[3:10] <- newnames[row.names(res1)[3:10]]
res1 <- round(res1,2)

cl2 <- newdat[kc$cluster==2,]
res2 <- data.frame(M=unlist(lapply(cl2[,vars],function(c) mean(c,na.rm=T))))
res2$CI_L <- sapply(vars
                    ,function(c) confint(lm(data=cl2,as.formula(paste0(c,"~1"))))[1] )
res2$CI_U <- sapply(vars
                    ,function(c) confint(lm(data=cl2,as.formula(paste0(c,"~1"))))[2] )
res2$SD = unlist(lapply(cl2[,vars],function(c) sd(c,na.rm=T)))
row.names(res2)[3:10] <- newnames[row.names(res2)[3:10]]
res2 <- round(res2,2)

df = data.frame(Cluster1=
apply(res1,1,function(x) paste0("M = ",x[1]," (",x[2],"-",x[3],")\nSD = ",x[4]))
,Cluster2=
  apply(res2,1,function(x) paste0("M = ",x[1]," (",x[2],"-",x[3],")\nSD = ",x[4])) )


tmp <- rbind(table(kc$cluster),table(newdat$casual,kc$cluster))
tmp <- as.data.frame(tmp)
names(tmp) <- c("Cluster1","Cluster2")
df <- rbind(tmp,df)

write.csv(df,"CompareClusters.csv")



# How other domains predict TPACK
f <- as.formula(paste0("P5~",paste(c(paste0("P",c(1:4,6:8)),"age_n","casual"),collapse="+")))
fitTPACK <- lm(f,data=newdat)
cf = summary(fitTPACK)$coef
row.names(cf)[2:8] <- newnames[row.names(cf)[2:8]]

write.csv(round(cf,3),"TPACK.csv")

newdat$confclust <- kc$cluster==1
fit_logit=glm(confclust~age_n+casual+expy,data=newdat,family=binomial)
summary(fit_logit)


fint <- as.formula(paste0("P5~",paste(paste(c(paste0("P",c(1:4,6:8)),"age_n","casual"),"confclust",sep="*"),collapse="+")))
fitTPACKint <- lm(fint,data=newdat)
summary(fitTPACKint)


library(relaimpo)

fitTPACK_conf <- lm(f,data=newdat,subset=confclust)
summary(fitTPACK_conf)
library(QuantPsyc)
beta=lm.beta(fitTPACK_conf)
# Best predictor of high TPACK is not-casual,low PK,high CK, high BNCL
ri_conf <- calc.relimp(fitTPACK_conf)$lmg
out=100*ri_conf[order(ri_conf,decreasing = T)]
df <- data.frame(predictors=names(out),RI=NA,B=NA)
tmp <- newnames[names(out)]
cf <- fitTPACK_conf$coef
names(cf) <- gsub("TRUE","",names(cf))
df <- data.frame(RI=round(out,1),B=round(cf[match(names(out),names(cf))],2) )
row.names(df)=ifelse(is.na(tmp),names(out),tmp)
write.csv(df,"conf.csv")


cf2 <- as.data.frame( summary(fitTPACK_conf)$coef[,c(1,2,4)] )
cf2 <- round(cf2,2)
cf2$beta <- c("",round(lm.beta(fitTPACK_conf),2))
cf2$relimp <- c("",round(ri_conf*100,1))
cf2 <- cf2[c(1,order(ri_conf,decreasing=T)+1),]
row.names(cf2) <- ifelse(is.na(newnames[row.names(cf2)]),row.names(cf2),newnames[row.names(cf2)])
write.csv(cf2,"cluster1_TPACK_reg.csv")

fitTPACK_notconf <- lm(f,data=newdat,subset=!confclust)
summary(fitTPACK_notconf)
# Best predictor of high TPACK is high TPK, low casual, high CK, low TCK
ri_nconf <- calc.relimp(fitTPACK_notconf)$lmg
out=100*ri_nconf[order(ri_nconf,decreasing = T)]
df <- data.frame(predictors=names(out),RI=NA,B=NA)
tmp <- newnames[names(out)]
cf <- fitTPACK_notconf$coef
names(cf) <- gsub("TRUE","",names(cf))
df <- data.frame(RI=round(out,1),B=round(cf[match(names(out),names(cf))],2) )
row.names(df)=ifelse(is.na(tmp),names(out),tmp)
write.csv(df,"nconf.csv")



cf2 <- as.data.frame( summary(fitTPACK_notconf)$coef[,c(1,2,4)] )
cf2 <- round(cf2,2)
cf2$beta <- c("",round(lm.beta(fitTPACK_notconf),2))
cf2$relimp <- c("",round(ri_nconf*100,1))
cf2 <- cf2[c(1,order(ri_nconf,decreasing=T)+1),]
row.names(cf2) <- ifelse(is.na(newnames[row.names(cf2)]),row.names(cf2),newnames[row.names(cf2)])
write.csv(cf2,"cluster2_TPACK_reg.csv")






png("CKage_TPACKexp.png",width=480,height=240)
par(mfrow=c(1,2))
plot(P6~age,data=newdat,ylab="CK",xlab="age",las=2)
plot(P5~expy,data=newdat,ylab="TPACK",xlab="years of experience",las=2)
dev.off()


write.csv(table(kc$cluster,newdat$age),"cluster_v_age.csv")




### CK Linear regression.
f <- as.formula(paste0("P6~",paste(c(paste0("P",c(1:5,7:8)),"age_n","casual"),collapse="+")))

fitCK_conf <- lm(f,data=newdat,subset=confclust)
summary(fitCK_conf)
beta=lm.beta(fitCK_conf)
ri_conf <- calc.relimp(fitCK_conf)$lmg

cf2 <- as.data.frame( summary(fitCK_conf)$coef[,c(1,2,4)] )
cf2 <- round(cf2,2)
cf2$beta <- c("",round(lm.beta(fitCK_conf),2))
cf2$relimp <- c("",round(ri_conf*100,1))
cf2 <- cf2[c(1,order(ri_conf,decreasing=T)+1),]
row.names(cf2) <- ifelse(is.na(newnames[row.names(cf2)]),row.names(cf2),newnames[row.names(cf2)])
cf2_conf <- cf2

fitCK_notconf <- lm(f,data=newdat,subset=!confclust)
summary(fitCK_notconf)
# Best predictor of high TPACK is high TPK, low casual, high CK, low TCK
ri_nconf <- calc.relimp(fitCK_notconf)$lmg

cf2 <- as.data.frame( summary(fitCK_notconf)$coef[,c(1,2,4)] )
cf2 <- round(cf2,2)
cf2$beta <- c("",round(lm.beta(fitCK_notconf),2))
cf2$relimp <- c("",round(ri_nconf*100,1))
cf2 <- cf2[c(1,order(ri_nconf,decreasing=T)+1),]
row.names(cf2) <- ifelse(is.na(newnames[row.names(cf2)]),row.names(cf2),newnames[row.names(cf2)])

cfCK <- rbind(cf2_conf,cf2)

write.csv(cfCK,"CK_reg.csv")