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, Physical_Aggression,  Physical_victimization)

Find Best Fitting Model

compare_solutions(data, SelfHarm_total, Physical_Aggression,  Physical_victimization, 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, Physical_Aggression,  Physical_victimization,
                        variances = "equal",
                        covariances = "zero",
                        n_profiles = 3,
                        return_orig_df = TRUE,
                        print_which_stats = "all")
## Fit NA model with 3 profiles.
## LogLik is 1201.161
## AIC is 2430.322
## CAIC is 2495.652
## BIC is 2481.652
## SABIC is 2437.256
## ICL is 2535.396
## Entropy is 0.922

Best Fitting

model2 <- estimate_profiles(data,  
                         SelfHarm_total, Physical_Aggression,  Physical_victimization,
                        variances = "equal",
                        covariances = "equal",
                        n_profiles = 3,
                        return_orig_df = TRUE,
                        print_which_stats = "all")
## Fit NA model with 3 profiles.
## LogLik is 1194.463
## AIC is 2422.925
## CAIC is 2502.255
## BIC is 2485.255
## SABIC is 2431.345
## ICL is 2535.68
## Entropy is 0.928
write.csv(model2, "model2_selfHarm_physAggression_physicalVictim.csv", row.names = F)
 model2
## # A tibble: 289 x 5
##    SelfHarm_total Physical_Aggres… Physical_victim… profile posterior_prob
##             <dbl>            <dbl>            <dbl> <fct>            <dbl>
##  1              0             3.17             3.5  1                0.999
##  2              3             1.33             2.5  2                0.977
##  3              0             1                2.38 2                0.980
##  4              0             1.17             1.38 2                1.000
##  5              1             1.17             1.75 2                0.999
##  6              1             1                1.12 2                1.000
##  7              0             3.83             3    1                1    
##  8              0             2.33             3    2                0.562
##  9              0             3                3    1                0.996
## 10              0             1.17             1.5  2                1.000
## # ... with 279 more rows
# 1, 1, 0, 0, 1, -, 1
# 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