library(party)
## Loading required package: grid
## Loading required package: mvtnorm
## Loading required package: modeltools
## Loading required package: stats4
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: sandwich
library(pspearman)

write smi as excel document and add the weights at the last row.

write.xlsx(smi, file="smi.xlsx", sheetName="Sheet1")
t5 <- read.csv("smi.csv", header = TRUE)


t5[3:66]<- data.frame(lapply(t5[3:66], function(X) X/X[2292]))
t5 = t5[-2292,]
t5[,3:66]<-scale(t5[,3:66]) 
t5[,70:72]<-scale(t5[,70:72]) 

t5$b[t5$swl >= 1 & t5$swl < 2.8 ] <- "one"  
t5$b[t5$swl >=  2.8 & t5$swl < 4.2 ] <- "two"  
t5$b[t5$swl >= 4.2  & t5$swl <= 7] <- "three" 

ind = sample(2, nrow(t5), replace = TRUE, prob=c(0.7, 0.3))
trainset1 = t5[ind == 1,]
testset1 = t5[ind == 2,]

training model: Please repeat the training model for more than 10 times and test the criterion validity when ACC > 0.58

fit1 <- cforest((b == 'three')~   affect + negemo+ future+swear+sad
                +negate+sexual+death + filler+leisure + conj+ funct +discrep + i
                +future + past + bio + body + cogmech + death + cause + quant
                + future +incl + motion  + sad + tentat + excl+insight +percept +posemo
                +ppron +quant + relativ + space + article + age + s_all + s_sad + gender
                , data = trainset1, 
                controls=cforest_unbiased(ntree=500, mtry= 1))

table1 <- table(predict(fit1, OOB=TRUE, type = 'response') > 0.5, trainset1$b == 'three') 

# Training Step 2. 

nonthree <- trainset1[trainset1$b != 'three',]
nonthree$b <- nonthree$b == 'two'
fit2 <- cforest(b~ affect +bio+body + cogmech + discrep + excl
                + i + insight + negemo + article + relativ +space + swear + percept
                + posemo + ppron + ipron + preps + pronoun + motion 
                + tentat + we + you+ conj+ funct +age + shehe + they
                + quant + s_all + s_sad + gender
                + future+swear+sad + ingest
                +sexual+death + filler+leisure 
                +incl + motion + quant + sad + tentat
                + family + cause + see + hear + feel+ verb + past + present
                + humans
                ,  data = nonthree,
                controls=cforest_unbiased(ntree=500, mtry = 1))

### here I adjust the probablity, only p> 0.6 will be classified as a two  (two = 521, one = 121)
table2<-table(predict(fit2, OOB=TRUE, type = 'response') > 0.6, nonthree$b)

#training step 3, predict "one" and "three" after step two

fit3 <- cforest((b == 'one')~  affect + negemo+ future+swear+sad
                +negate+sexual+death + filler+leisure + conj+ funct +discrep + i
                +future + past + bio + body + cogmech + death + cause + quant
                +incl + motion  + sad + tentat + excl+insight +percept +posemo
                +ppron +quant + relativ + space + article + age + s_all + s_sad + gender 
                , data = trainset1, 
                controls=cforest_unbiased(ntree=500, mtry= 1))
# Testing with separate test set (testset2)

#run the test set through step1 
test_predict1<-predict(fit1, newdata=testset1, type='response')
#produce the table for splitting 'three' and 'other'
test_table1<-table(test_predict1>0.5, testset1$b=="three")

#get the data that is classed as 'other', this will be passed to step 2
step2_index<-test_predict1<=0.5

#create the dataset that would be passed to step 2 (this will probably contain 'one', 'two' and 'three' data)
testset2<-testset1[step2_index,]
# run the second step on test set 2
test_predict2<-predict(fit2, newdata=testset2, type='response')
#find out how many 'two' would be correctly classified
test_table2a<-table(test_predict2>0.6, testset2$b=="two")

step3_index<-test_predict2<=0.6
testset3<-testset2[step3_index,]

test_predict3<-predict(fit3, newdata=testset3, type='response')
#find out how many 'one' would be correctly classified
#here I justed the threshold to 0.23, because the response in this step ranged from 0.1-0.3
test_table3a<-table(test_predict3 <= 0.23, testset3$b=="three")


# the equation for the test set accuracy is then:
total_percentage_test<-(test_table1[2,2]+test_table2a[2,2]+test_table3a[1,1]+test_table3a[2,2])/nrow(testset1)*100 

total_percentage_train<-(table1[2,2]+table2[2,2]+table2[1,1])/nrow(trainset1)*100
# disclosure

smi <- read.csv("smi.csv", header = TRUE)
sdfm <- read.csv("sdfm.csv", header = TRUE)
s_merged <- merge(sdfm, smi, by.x="userid", by.y="userid", all.x=TRUE)
t14 <- s_merged[complete.cases(s_merged),]

#write.xlsx(t14, file="t14.xlsx", sheetName="Sheet1")

t5 <- read.csv("t14.csv", header = TRUE)


t5 <- as.data.frame(t5)

t5[13:76]<- data.frame(lapply(t5[13:76], function(X) X/X[896]))
t5 = t5[-896,]
t5[,13:76]<-scale(t5[,13:76]) 
t5[,70:72]<-scale(t5[,70:72]) 

t5$b[t5$swl >= 1 & t5$swl < 2.8 ] <- "1"  
t5$b[t5$swl >=  2.8 & t5$swl < 4.2 ] <- "2"  
t5$b[t5$swl >= 4.2  & t5$swl <= 7] <- "3" 
t5$b <- as.numeric(t5$b)  

testset2 <- t5

cor1 <- cor.test(testset2$sd_score, testset2$b, method="spearman", exact=FALSE)
cor2 <- cor.test(testset2$fm_score, testset2$b, method="spearman", exact=FALSE)
cor3 <- cor.test(testset2$work.x, testset2$b, method="spearman", exact=FALSE)

cor1
## 
##  Spearman's rank correlation rho
## 
## data:  testset2$sd_score and testset2$b
## S = 102270000, p-value = 1.514e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1440697
cor2
## 
##  Spearman's rank correlation rho
## 
## data:  testset2$fm_score and testset2$b
## S = 105340000, p-value = 0.0003873
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1183647
cor3
## 
##  Spearman's rank correlation rho
## 
## data:  testset2$work.x and testset2$b
## S = 121770000, p-value = 0.5685
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.01908596
testset2$pre_swl<-predict(fit1, newdata=testset2 , type='response')
testset2$pre_swl[testset2$pre_swl > 0.5] <- 3

#get the data that is classed as 'other', this will be passed to step 2
step2_index<-testset2$pre_swl<= 0.5
testset3<-testset2[step2_index,]

testset3$pre_swl3<-predict(fit2, newdata=testset3, type='response')
testset3$pre_swl3[testset3$pre_swl3> 0.6] <- 2

step3_index<-test_predict2<= 0.6
testset4<-testset3[step3_index,]

#Step three 
testset4$pre_swl4<-predict(fit3, newdata=testset4, type='response')
testset4$pre_swl4[testset4$pre_swl4 <= 0.23] <- 3
testset4$pre_swl4[testset4$pre_swl4 > 0.23 & testset4$pre_swl4 < 3 ] <- 1

testset5 <- merge(testset2, testset3, by = "userid", all.x=T)
testset5$swl2 <- ifelse(is.na(testset5$pre_swl3), testset5$pre_swl.x, testset5$pre_swl3)
testset6 <- merge(testset5, testset4, by = "userid", all.x=T)
testset6$swl3 <- ifelse(is.na(testset6$pre_swl4), testset6$swl2, testset6$pre_swl4)
testset6$swl3[testset6$swl3 < 1 ] <- 2


cor1 <- cor.test(testset6$sd_score.x, testset6$swl3, method="spearman", exact=FALSE)
cor2 <- cor.test(testset6$fm_score.x, testset6$swl3, method="spearman", exact=FALSE)
cor3 <- cor.test(testset6$work.x, testset6$swl3, method="spearman", exact=FALSE)

cor1
## 
##  Spearman's rank correlation rho
## 
## data:  testset6$sd_score.x and testset6$swl3
## S = 178220000, p-value = 6.122e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1855644
cor2
## 
##  Spearman's rank correlation rho
## 
## data:  testset6$fm_score.x and testset6$swl3
## S = 258110000, p-value = 2.199e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.1795221
cor3
## 
##  Spearman's rank correlation rho
## 
## data:  testset6$work.x and testset6$swl3
## S = 1782700, p-value = 1.522e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3843446