Load Libraries

library(foreign)
library(tidyLPA)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
options(scipen = 999)

Load Data

data <- read.spss("data_306.sav", to.data.frame = T)
## Warning in read.spss("data_306.sav", to.data.frame = T): data_306.sav:
## Very long string record(s) found (record type 7, subtype 14), each will be
## imported in consecutive separate variables
## Warning in read.spss("data_306.sav", to.data.frame = T): Undeclared
## level(s) 2 added in variable: profile
## Warning in read.spss("data_306.sav", to.data.frame = T): Undeclared
## level(s) 0 added in variable: profile_dummy
#data <- data %>% filter(., Gender == "אישה")   #
        
#"גבר"

#"אישה"

# data <- select(data, serial, Abuse, Gender, RTA_factor_1, RTA_factor_2,
#                RTA_factor_3, RTA_factor_4, SexVictimization_no1,
#                Physical_victimization, SexVictimization_agg_no1,
#                Physical_Aggression, Dissociation_total, SelfHarm_total, guilt_tot,
#                SexVictimization_agg_no1, Emotional_abuse)

data <- filter(data, Abuse == "reported an abuse")

#data <- select(data, RTA_factor_1, RTA_factor_2, RTA_factor_3, RTA_factor_4)
data <- select(data, SelfHarm_total, SexVictimization_no1, SexVictimization_agg_no1)

Find Best Fitting Model

compare_solutions(data, SelfHarm_total, SexVictimization_no1, SexVictimization_agg_no1,
                  scale_raw_data = T)

It appears as though model 2 with 2 profiles fits best (as indicated by low BIC). Model 3 with 3 profiles has the lowest BIC but I’m suspicious of it due to large variance.

Look at all stats

model1 <- estimate_profiles(data,  
                        SelfHarm_total, SexVictimization_no1, SexVictimization_agg_no1,
                        variances = "equal",
                        covariances = "zero",
                        n_profiles = 3,
                        return_orig_df = TRUE,
                        print_which_stats = "all")
## Fit NA model with 3 profiles.
## LogLik is 1549.229
## AIC is 3126.458
## CAIC is 3190.361
## BIC is 3176.361
## SABIC is 3131.975
## ICL is 3193.283
## Entropy is 0.972

Best Fitting

model2 <- estimate_profiles(data,  
                         SelfHarm_total, SexVictimization_no1, SexVictimization_agg_no1,
                        variances = "equal",
                        covariances = "equal",
                        n_profiles = 2,
                        return_orig_df = TRUE,
                        print_which_stats = "all")
## Fit NA model with 2 profiles.
## LogLik is 1593.449
## AIC is 3212.899
## CAIC is 3272.238
## BIC is 3259.238
## SABIC is 3218.022
## ICL is 3273.521
## Entropy is 0.975
write.csv(model2, "model2.csv", row.names = F)
model2
## # A tibble: 261 x 5
##    SelfHarm_total SexVictimizatio… SexVictimizatio… profile posterior_prob
##             <dbl>            <dbl>            <dbl> <fct>            <dbl>
##  1              0                0                0 1                1.000
##  2              3                4                0 1                0.904
##  3              1                2                3 1                1.000
##  4              1                2                0 1                1.000
##  5              0                0                0 1                1.000
##  6              0                4                2 1                1.000
##  7              0                6                6 1                1.000
##  8              0                1                1 1                1.000
##  9              1                5                2 1                1.000
## 10              0                2                1 1                1.000
## # ... with 251 more rows
# model3 <- estimate_profiles(data,  
#                         SelfHarm_total, Physical_Aggression, Physical_victimization,
#                         variances = "varying",
#                         covariances = "zero",
#                         n_profiles = 2,
#                         return_orig_df = TRUE,
#                         print_which_stats = "all")
# model6 <- estimate_profiles(data,  
#                         SelfHarm_total, Physical_Aggression, Physical_victimization,
#                         variances = "varying",
#                         covariances = "varying",
#                         n_profiles = 2,
#                         return_orig_df = TRUE,
#                         print_which_stats = "all")
# bootstrap_lrt(data, RTA_factor_1, RTA_factor_2, 
#               RTA_factor_3, RTA_factor_4, model = 1)

y-axis = mean

plot_profiles(model2) #to_center = TRUE, to_scale = TRUE

y-axis = mean centered

print(plot_profiles(model2, to_center = T)) #to_center = TRUE, to_scale = TRUE

y-axis = standard deviation

plot_profiles(model2, to_center = T, to_scale = T) #to_center = TRUE, to_scale = TRUE