## data read
daf <- read.csv('D:/all_variables2.csv')

head(daf)
dim(daf)

#daf <- subset(daf,TC95>0.95)
summary(daf)
das <- daf[sample(nrow(daf),600),]
cor.test(das$FSDI,das$SIF)
summary(das)

hist(das$FSDI)
hist(das$Richness)
#das$FSDI <- log(das$FSDI)

#fit <-  lm(Richness~ RAD+PRE + TMP +POP + FIRE+ DEM+CEC + SOC + SAND , das)
#vif(fit)


rich <- lm(Richness~ RAD+PRE + TMP +POP + FIRE+ DEM+CEC + SOC + SAND, das )

fsdi <- lm(FSDI ~ RAD+PRE + TMP +POP + FIRE+ DEM+CEC + SOC + SAND,das )

gpp <- lm(SIF~FSDI + Richness +RAD+PRE + TMP +POP + FIRE+ DEM+CEC + SOC + SAND ,das)


model <- psem(rich, fsdi, gpp)

fisherC(model)

summary(model)

plot(model)


anova(model)

write.csv(das, "D:/ALL_stats_TC95_sif600point.csv")