#######################   load R packages   ####################################

library("metaSEM")
library("foreign")
library("semPlot")
library("lavaan")

#######################   Clean the data   ####################################

dataset=read.spss("01 TSSEM_komplett.sav", use.value.labels=FALSE, to.data.frame=T)

head(dataset)

##Leere Felder löschen
dataset <- dataset[-c(177, 178, 179, 180, 181, 182, 183, 184), ]
dataset



##Missings setzen
dataset$Zelle_1[dataset$Zelle_1==-88] <- NA
dataset$Zelle_2[dataset$Zelle_2==-88] <- NA
dataset$Zelle_3[dataset$Zelle_3==-88] <- NA
dataset$Zelle_4[dataset$Zelle_4==-88] <- NA
dataset$Zelle_5[dataset$Zelle_5==-88] <- NA
dataset$Zelle_6[dataset$Zelle_6==-88] <- NA
dataset$Zelle_7[dataset$Zelle_7==-88] <- NA
dataset$Zelle_8[dataset$Zelle_8==-88] <- NA
dataset$Zelle_9[dataset$Zelle_9==-88] <- NA
dataset$Zelle_10[dataset$Zelle_10==-88] <- NA
dataset$Zelle_1[dataset$Zelle_1==-99] <- NA
dataset$Zelle_2[dataset$Zelle_2==-99] <- NA
dataset$Zelle_3[dataset$Zelle_3==-99] <- NA
dataset$Zelle_4[dataset$Zelle_4==-99] <- NA
dataset$Zelle_5[dataset$Zelle_5==-99] <- NA
dataset$Zelle_6[dataset$Zelle_6==-99] <- NA
dataset$Zelle_7[dataset$Zelle_7==-99] <- NA
dataset$Zelle_8[dataset$Zelle_8==-99] <- NA
dataset$Zelle_9[dataset$Zelle_9==-99] <- NA
dataset$Zelle_10[dataset$Zelle_10==-99] <- NA
dataset$WR[dataset$WR==-999] <- NA



#######################  Correlation matrices     #######################

nvar <- 5
varnames <- c("Religiosität", "soziale", "Optimismus", "Selbstwert", "Lebenszufriedenheit")
labels <- list (varnames, varnames)


cormatrices <- list ()

for (i in 1:nrow(dataset)) {
  cormatrices [[i]]<-vec2symMat(as.matrix(dataset[i,2:11]), 
                                diag = FALSE)
  dimnames(cormatrices[[i]])<-labels
}

head(cormatrices)


#######################   OSMASEM ####################################



#Cor matrices in Data-Frame 

my.df <- Cor2DataFrame(cormatrices, dataset$Fallzahl)
my.df

## Add the standardized HDI (moderator) to the data
my.df$data <- data.frame(my.df$data, HDI_Stand=scale(dataset$HDI), check.names=FALSE)
head(my.df$data)

## Display the pairwise no. of studies
pattern.na(my.df, show.na=FALSE, type="osmasem")

#Modell in lavaan-Code

vollModell <-'soziale ~ b31*Religiosität
Optimismus ~ b32*Religiosität
Selbstwert ~ b33*Religiosität
Lebenszufriedenheit ~ b34*Religiosität
Lebenszufriedenheit ~ b35*soziale
Lebenszufriedenheit ~ b36*Optimismus
Lebenszufriedenheit ~ b37*Selbstwert
soziale ~~ k1*Optimismus
soziale ~~ k2*Selbstwert
Optimismus ~~ k3*Selbstwert
Religiosität ~~ 1*Religiosität
soziale ~~ p22*soziale
Optimismus ~~ p33*Optimismus
Selbstwert ~~ p44*Selbstwert
Lebenszufriedenheit ~~p55*Lebenszufriedenheit

'


plot(vollModell, col="grey")    

#Transform into RAM
RAM1 <- lavaan2RAM(vollModell, obs.variables=c("Religiosität", "soziale", "Optimismus", "Selbstwert", "Lebenszufriedenheit"))
RAM1



## Create the model implied correlation structure with implicit diagonal constraints
M0 <- create.vechsR(A0=RAM1$A, S0=RAM1$S)

## Create the heterogeneity variance-covariance matrix
T0 <- create.Tau2(RAM=RAM1, RE.type="Diag")

## Fit the OSMASEM
fit0 <- osmasem(model.name="Just identified model", Mmatrix=M0, Tmatrix=T0, data=my.df)
summary(fit0)

## Get the estimated coefficients
coef(fit0)

## A matrix
mxEval(Amatrix, fit0$mx.fit)

## S matrix
mxEval(Smatrix, fit0$mx.fit)

## Extract the heterogeneity variance-covariance matrix
VarCorr(fit0)

#A Matrix alle Regressionskoeffizienten
Ax<-create.mxMatrix(c
                    (0, 0, 0, 0, 0,
                      "0*data.HDI", 0, 0, 0, 0,
                      "0*data.HDI", 0, 0, 0, 0,
                      "0*data.HDI", 0, 0, 0, 0,
                      "0*data.HDI", "0*data.HDI", "0*data.HDI", "0*data.HDI", 0), 
                    type="Full",
                    ncol=5,
                    nrow=5,
                    byrow = TRUE)
Ax

## Create the model implied correlation structure with the centered HDI as the moderator

M1 <- metaSEM::create.vechsR(A0=RAM1$A, S0=RAM1$S, Ax=Ax)

fit1 <- osmasem(model.name="HDI as moderator on Ax", Mmatrix=M1, Tmatrix=T0, data=my.df)
summary(fit1)


#######################   Ende der Analyse ####################################
sessionInfo(package = NULL)


######################################################################################################################################



my.df <- Cor2DataFrame(Roorda11$data, Roorda11$n)

## Add centered SES as a moderator (standardizing SES helps to improve the stability of the results)
my.df$data <- data.frame(my.df$data, SES=scale(Roorda11$SES), check.names=FALSE)
head(my.df$data)



pattern.na(my.df, show.na=FALSE, type="osmasem")



## Proposed model
model2 <- 'enga ~ b31*neg + b32*pos
           achiev ~ b41*neg + b42*pos + b43*enga
           pos ~~ p21*neg
           pos ~~ 1*pos
           neg ~~ 1*neg
           enga ~~ p33*enga
           achiev ~~ p44*achiev'

plot(model2, layout="spring", color="yellow")

## Convert the lavaan model to the RAM specification
RAM2 <- lavaan2RAM(model2, obs.variables=c("neg", "pos", "enga", "achiev"))
RAM2

## Create matrices with implicit diagonal constraints
M0 <- create.vechsR(A0=RAM2$A, S0=RAM2$S)

## Create heterogeneity variances
T0 <- create.Tau2(RAM=RAM2, RE.type="Diag")

mx.fit0a <- osmasem(model.name="Just identified model", Mmatrix=M0, Tmatrix=T0, data= my.df)
summary(mx.fit0a)

coef(mx.fit0a)


## Extract the variance component
VarCorr(mx.fit0a)


## Proposed model
model3 <- 'enga ~ b31*neg + b32*pos
           achiev ~ 0*neg + 0*pos + b43*enga
           pos ~~ p21*neg
           pos ~~ 1*pos
           neg ~~ 1*neg
           enga ~~ p33*enga
           achiev ~~ p44*achiev'

plot(model3, layout="spring", color="yellow")



RAM3 <- lavaan2RAM(model3, obs.variables=c("neg", "pos", "enga", "achiev"))
RAM3

## Create matrices with implicit diagonal constraints
M0b <- create.vechsR(A0=RAM3$A, S0=RAM3$S)

## Create heterogeneity variances
T0b <- create.Tau2(RAM=RAM3, RE.type="Diag")

mx.fit0b <- osmasem(model.name="Over identified model", Mmatrix=M0b, Tmatrix=T0b, data= my.df)

## Get the summary and test statistics
summary(mx.fit0b, Saturated=TRUE)

## Get the SRMR
osmasemSRMR(mx.fit0b)

## Get the coefficients
coef(mx.fit0b)

## Extract the variance component
VarCorr(mx.fit0b)

## data.SES as a moderator
Ax <- matrix(c(0,0,0,0,
               0,0,0,0,
               "0*data.SES","0*data.SES",0,0,
               "0*data.SES","0*data.SES","0*data.SES",0),
             nrow=4, ncol=4, byrow=TRUE)
Ax

M1 <- create.vechsR(A0=RAM2$A, S0=RAM2$S)

mx.fit1 <- osmasem(model.name="Ax as moderator", Mmatrix=M1, Tmatrix=T0, data= my.df)
summary(mx.fit1)

## Get the coefficients
coef(mx.fit1)



## Get the R2
osmasemR2(mx.fit1, mx.fit0a)


## Extract the variance component
VarCorr(mx.fit1)



## Comparing with the model without the moderator
anova(mx.fit1, mx.fit0a)