LVTU PSC-17 analüüsid

Andmed

# some simpler variable names
D$att <- D$`Y_PSC-17.Attention.Subscale`
D$int <- D$`Y_PSC-17.Internalizing.Subscale`
D$ext <- D$`Y_PSC-17.Externalizing.Subscale`
D$attlv <- D$`LV_PSC-17.Attention.Subscale`
D$intlv <- D$`LV_PSC-17.Internalizing.Subscale`
D$extlv <- D$`LV_PSC-17.Externalizing.Subscale`
D$sugu2 <- factor(D$sugu, levels=1:2)
D$vg2 <- factor(D$vanus_grupp, levels=2:3)

psc_laps <- c("Y_PSC_1", "Y_PSC_2", "Y_PSC_3", "Y_PSC_4", "Y_PSC_5", "Y_PSC_6",  "Y_PSC_7", "Y_PSC_8", "Y_PSC_9", "Y_PSC_10", "Y_PSC_11", "Y_PSC_12", "Y_PSC_13", "Y_PSC_14", "Y_PSC_15", "Y_PSC_16", "Y_PSC_17")
psc_vanem <- gsub("Y", "LV", psc_laps)
names(D)[names(D) %in% psc_laps] <- paste0("PSC", 1:17) -> psc_l
names(D)[names(D) %in% psc_vanem] <- paste0("LV", 1:17) -> psc_lv
att <- c(1,3,7,13,17)
int <- c(2,6,9,11,15)
ext <- setdiff(1:17, c(att,int))

with(D, fftable(ankeedikeel, LV_ankeedikeel, exclude=NULL, labs=c("Laps", "Lapsevanem")))

	Lapsevanem
Laps	1	2	Kokku
1	533	10	543
2	15	123	138
Kokku	548	133	681

incl <- D$ankeedikeel == 1 | D$LV_ankeedikeel == 1
Di <- subset(D, ankeedikeel == 1)
Dj <- subset(D, LV_ankeedikeel== 1)
Dij <- subset(D, ankeedikeel == 1 & LV_ankeedikeel==1)
D <- subset(D, incl)

Jätame andmestikku vastajad, kelle puhul kas lapsevanem või laps on vastanud eesti keeles. Sisse jäävad mõned üksikud 18- ja 19-aastased.

Valimi kirjeldus

Lapse sugu ja vanus (siin ainult lapsed, kes vastasid ankeedile eesti keeles)

fftable(Di$sugu, Di$vanus, labs = c("sugu", "vanus"))

	vanus
sugu	8	9	10	11	12	13	14	15	16	17	18	19	Kokku
1	13	33	33	25	26	27	20	26	24	25	6	0	258
2	13	36	24	29	20	24	23	28	38	33	10	1	279
3	0	1	0	0	0	0	0	0	0	3	0	0	4
4	0	0	0	0	1	0	0	0	0	1	0	0	2
Kokku	26	70	57	54	47	51	43	54	62	62	16	1	543

Lapse sugu ja vanuserühm (1= 8-10; 2= 11-14; 3= 15+)

fftable(Di$sugu, Di$vanus_grupp, labs = c("sugu", "vanuserühm"))

	vanuserühm
sugu	1	2	3	Kokku
1	79	98	81	258
2	73	96	110	279
3	1	0	3	4
4	0	1	1	2
Kokku	153	195	195	543

Lapsevanema sugu ja vanus (ainult lapsevanemad, kes vastasid ankeedile eesti keeles)

Dj$LV_vg <- cut(Dj$LV_vanus, 2:7*10)
fftable(Dj$LV_sugu, Dj$LV_vg, labs=c("sugu", "vanus"))

	vanus
sugu	(20,30]	(30,40]	(40,50]	(50,60]	(60,70]	Kokku
1	0	10	21	5	1	37
2	9	154	184	36	1	384
Kokku	9	164	205	41	2	421

• Laps & lapsevanem vastasid mõlemad eesti keeles: N = 533
• Kas laps või lapsevanem vastas eesti keeles: N = 558

# Korrelatsioonid 2. vs 3. vanuserühmas
c1 <- cor(subset(Di, Di$vanus_grupp==2, select=psc_l), use="pai")
c2 <- cor(subset(Di, Di$vanus_grupp==3, select=psc_l), use="pai")
plot(lt(c1), lt(c2), xlab="Korrelatsioonid 2. vanuserühmas", ylab="Korrelatsioonid 3. vanuserühmas")

CFA models and invariance testing

Choices in CFA

Lühike lõik, mida saab kasutada meetodi-osas:

We used the WLSMV estimator due to its suitability for ordered categorical indicators and its less stringent sample size requirements compared to alternatives (Brown, 2015). Following the recommendations of Brosseau-Liard and Savalei (2014), we applied robust corrections to the fit indices (CFI, TLI, RMSEA, and SRMR). To interpret the invariance tests, we used the following cutoff values for noninvariance based on Chen (2007): (1) for metric invariance (testing for equality of loadings), a CFI decrease of .010 or greater along with an RMSEA increase of .015 or greater or an SRMR increase of .030 or greater; (2) for scalar invariance (testing for equality of intercepts), a CFI decrease of .010 or greater combined with an RMSEA increase of .015 or greater or an SRMR increase of .010 or greater. Scaled chi-square and df are reported.

Brosseau-Liard, P. E., & Savalei, V. (2014). Adjusting Incremental Fit Indices for Nonnormality. Multivariate Behavioral Research, 49(5), 460–470. https://doi.org/10.1080/00273171.2014.933697

Brown, T. A. (2015). Confirmatory factor analysis for applied research (2nd ed.). Guilford Press. [ ISBN 978-1-4625-1536-3]

Chen, F. F. (2007). Sensitivity of goodness of fit indexes to lack of measurement invariance. Structural Equation Modeling, 14(3), 464–504. https://doi.org/10.1080/10705510701301834

Quotes for justifying the method

• Some examples of invariance testing: https://dlab.berkeley.edu/news/testing-measurement-invariance-using-lavaan-r

1. Using WLSMV estimator BECAUSE: “The WLSMV is a robust estimator which does not assume normally distributed variables and provides the best option for modelling categorical or ordered data (Brown, 2006)”.

Brown, T. (2006). Confirmatory factor analysis for applied research. New York: Guildford.

From 2015 edition (ISBN 978-1-4625-1536-3) “There are several estimators that can be used with categorical indicators: for example, WLS, robust WLS (WLSMV), and unweighted least squares (ULS).”

“At this writing, the Mplus program appears to provide the best options for CFA modeling with categorical data. This is due in part to the availability of the WLSMV estimator and the ability to estimate models with ML via numerical integration (ML estimation is discussed later in this chapter). The WLSMV estimator provides WLS parameter estimates by using a diagonal weight matrix (W) and robust standard errors, and a [Data Issues in CFA 355] mean- and variance-adjusted c2 test statistic (Muthén & Muthén, 1998–2012). Unlike WLS, WLSMV does not require W to be positive definite because W is not inverted as part of the estimation procedure. In WLSMV, the number of elements in the diagonal W equals the number of sample correlations in S, but this matrix is not inverted during estimation. Nevertheless, WLSMV estimation is fostered when N is larger than the number of rows in W. In the computation of the c2 test statistic and standard errors, the full W is used but not inverted. Thus, although CFA models using categorical indicators necessitate larger samples than comparably sized models using continuous indicators, the sample size requirements of WLSMV are far less restrictive than those of WLS. For example, Muthén has conducted unpublished simulation studies and has found that Ns of 150 to 200 may be sufficient for medium-sized models (e.g., 10–15 indicators). Flora and Curran (2004) have confirmed these results by showing that WLSMV produces accurate test statistics, parameter estimates, and standard errors of CFA models under a variety of conditions (e.g., sample sizes ranging from 100 to 1,000; varying degrees of non-normality and model complexity). In addition, preliminary simulation research has shown that WLSMV performs well with samples as low as 200 for variables with floor or ceiling effects (although, as with continuous indicators, very skewed categorical indicators call for larger Ns). More studies are needed to more fully establish the performance of WLSMV under various sample sizes and other conditions (e.g., skewness, model complexity, size of indicator relationships, behavior of goodness-of-fit statistics). Because requisite sample size is closely tied to the specific model and data of a given study, general rules of thumb are of limited utility. Instead, researchers are encouraged to use the Monte Carlo routines available in the recent versions of most latent variable programs to determine the necessary N for their particular CFA investigation (see Chapter 10).”

2. Using robust correction for fit indexes as recommended in https://doi.org/10.1080/00273171.2012.715252 https://pubmed.ncbi.nlm.nih.gov/26735008/ (“We show that the virtually unknown correction due to Li and Bentler (2006) , which we label the sample-corrected robust RMSEA, is a consistent estimate of the population ML RMSEA yet drastically reduces bias due to nonnormality in small samples.”) https://pubmed.ncbi.nlm.nih.gov/26732359/ https://doi.org/10.1080/00273171.2014.933697

3. CFI change for invariance tests

For sample sizes with adequate power, equal group sizes, and mixed invariance (i.e., some loadings are higher and some lower in the first group), Chen (2007) also suggested a criterion of a -.01 change in CFI, paired with changes in RMSEA of .015 and SRMR of .030 (for metric invariance) or .015 (for scalar or residual invariance). Source: https://pmc.ncbi.nlm.nih.gov/articles/PMC5145197/ https://doi.org/10.1016/j.dr.2016.06.004

Chen 2007: When sample size is small (total N ≤ 300), sample sizes are unequal, and the pattern of noninvariance is uniform, the following cutoff criteria are suggested: for testing loading invariance, a change of ≤ −.005 in CFI, supplemented by a change of ≥ .010 in RMSEA or a change of ≥ .025 in SRMR would indicate noninvariance; for testing intercept or residual invariance, a change of ≥ −.005 in CFI, supplemented by a change of ≥ .010 in RMSEA or a change of ≥ .005 in SRMR would indicate noninvariance.

When sample size is adequate (total N > 300) and sample sizes are equal across the groups, particularly when lack of invariance is mixed, more stringent criteria are suggested. For testing loading invariance, a change of ≥ −.010 in CFI, supplemented by a change of ≥ .015 in RMSEA or a change of ≥ .030 in SRMR would indicate noninvariance; for testing intercept or residual invariance, a change of ≥ −.010 in CFI, supplemented by a change of ≥ .015 in RMSEA or a change of ≥ .010 in SRMR would indicate noninvariance. Among the three indexes, CFI was chosen as the main criterion because RMSEA and SRMR tend to overreject an invariant model when sample size is small, particularly when using SRMR for testing loading or residual variance invariance. In addition, changes in RMSEA are more likely to be affected by sample size and model complexity. However, these criteria should be used with caution, because testing measurement invariance is a very complex issue.

Single-group model

mud1 <- " ATT =~ PSC1 + PSC3 + PSC7 + PSC13 + PSC17
         INT =~ PSC2 + PSC6 + PSC9 + PSC11 + PSC15
         EXT =~ PSC4 + PSC5 + PSC8 + PSC10 + PSC12 + PSC14 + PSC16
"
fit1 <- cfa(mud1, data=Di, estimator="WLSMV")
#fit12 <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2), estimator="WLSMV")
#fit13 <- cfa(mud1, data=subset(Di, vanus_grupp %in% 3), estimator="WLSMV")

Fit of single-group model in overall sample:

getfms(fit1)

chisq.sc	df.sc	pvalue.sc	cfi.r	tli.r	srmr	rmsea.r	rmsea.ci.lo.r	rmsea.ci.up.r
208	116	0.00000038	0.994	0.993	0.0512	0.0406	0.0315	0.0495

semPaths(fit1, rotation=2, what="std", residuals=FALSE, edge.color="black", esize=1, edge.label.cex=1, sizeMan = 6, node.width=1, node.height=1, weighted=FALSE, nodeLabels = c(PSC1 = "i1", PSC3 = "i3", PSC7 = "i7", PSC13 = "i13", PSC17 = "i17",  PSC2 = "i2", PSC6 = "i6", PSC9 = "i9", PSC11 = "i11", PSC15 = "i15", PSC4 = "i4", PSC5 = "i5", PSC8 = "i8", PSC10 = "i10", PSC12 = "i12", PSC14 = "i14", PSC16 = "i16", ATT = "ATT", INT = "INT", EXT = "EXT" ))

ss <- standardizedsolution(fit1) |> as.data.frame()
ss[,] <- lapply(ss, formatif, 3)
qflextable(ss)

lhs	op	rhs	est.std	se	z	pvalue	ci.lower	ci.upper
ATT	=~	PSC1	0.690	0.0360	19.19	0.0000000000000000	0.620	0.761
ATT	=~	PSC3	0.518	0.0419	12.38	0.0000000000000000	0.436	0.600
ATT	=~	PSC7	0.751	0.0277	27.12	0.0000000000000000	0.697	0.805
ATT	=~	PSC13	0.493	0.0422	11.67	0.0000000000000000	0.410	0.575
ATT	=~	PSC17	0.686	0.0342	20.07	0.0000000000000000	0.619	0.753
INT	=~	PSC2	0.736	0.0339	21.73	0.0000000000000000	0.670	0.803
INT	=~	PSC6	0.769	0.0290	26.51	0.0000000000000000	0.712	0.826
INT	=~	PSC9	0.646	0.0391	16.51	0.0000000000000000	0.569	0.723
INT	=~	PSC11	0.615	0.0425	14.46	0.0000000000000000	0.532	0.698
INT	=~	PSC15	0.678	0.0356	19.02	0.0000000000000000	0.608	0.747
EXT	=~	PSC4	0.486	0.0548	8.86	0.0000000000000000	0.378	0.593
EXT	=~	PSC5	0.548	0.0567	9.65	0.0000000000000000	0.437	0.659
EXT	=~	PSC8	0.323	0.0607	5.31	0.0000001070205649	0.204	0.442
EXT	=~	PSC10	0.499	0.0555	8.99	0.0000000000000000	0.390	0.608
EXT	=~	PSC12	0.537	0.0537	9.99	0.0000000000000000	0.431	0.642
EXT	=~	PSC14	0.389	0.0628	6.19	0.0000000005835605	0.266	0.512
EXT	=~	PSC16	0.357	0.0639	5.59	0.0000000230409289	0.232	0.482
PSC1	~~	PSC1	0.523	0.0497	10.53	0.0000000000000000	0.426	0.621
PSC3	~~	PSC3	0.731	0.0434	16.85	0.0000000000000000	0.646	0.816
PSC7	~~	PSC7	0.436	0.0416	10.47	0.0000000000000000	0.354	0.517
PSC13	~~	PSC13	0.757	0.0416	18.21	0.0000000000000000	0.676	0.839
PSC17	~~	PSC17	0.529	0.0469	11.29	0.0000000000000000	0.438	0.621
PSC2	~~	PSC2	0.458	0.0499	9.17	0.0000000000000000	0.360	0.556
PSC6	~~	PSC6	0.409	0.0446	9.16	0.0000000000000000	0.321	0.496
PSC9	~~	PSC9	0.583	0.0505	11.53	0.0000000000000000	0.484	0.682
PSC11	~~	PSC11	0.622	0.0523	11.89	0.0000000000000000	0.519	0.724
PSC15	~~	PSC15	0.541	0.0483	11.20	0.0000000000000000	0.446	0.636
PSC4	~~	PSC4	0.764	0.0532	14.36	0.0000000000000000	0.660	0.868
PSC5	~~	PSC5	0.700	0.0622	11.26	0.0000000000000000	0.578	0.822
PSC8	~~	PSC8	0.896	0.0392	22.86	0.0000000000000000	0.819	0.973
PSC10	~~	PSC10	0.751	0.0555	13.54	0.0000000000000000	0.642	0.859
PSC12	~~	PSC12	0.712	0.0577	12.35	0.0000000000000000	0.599	0.825
PSC14	~~	PSC14	0.849	0.0488	17.40	0.0000000000000000	0.753	0.945
PSC16	~~	PSC16	0.872	0.0457	19.11	0.0000000000000000	0.783	0.962
ATT	~~	ATT	1.000	0.0000	NA	NA	1.000	1.000
INT	~~	INT	1.000	0.0000	NA	NA	1.000	1.000
EXT	~~	EXT	1.000	0.0000	NA	NA	1.000	1.000
ATT	~~	INT	0.770	0.0358	21.49	0.0000000000000000	0.700	0.840
ATT	~~	EXT	0.680	0.0511	13.29	0.0000000000000000	0.579	0.780
INT	~~	EXT	0.462	0.0613	7.54	0.0000000000000464	0.342	0.583

Age group invariance

fit1mg <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2:3), estimator = "WLSMV", group="vanus_grupp")
sf1mg <- summary(fit1mg, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg)

#### Age group metric invariance (multigroup model, equal loadings)
fit1mg.el <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2:3), estimator = "WLSMV", group="vanus_grupp", group.equal="loadings")
sf1mg.el <- summary(fit1mg.el, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg.el)
cf1  <- compareFit(fit1mg, fit1mg.el)
t1 <- cfxtract(cf1, "metric invariance") 

#### Age group scalar invariance (multigroup model, equal intercepts) 
fit1mg.ei <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2:3), estimator = "WLSMV", group="vanus_grupp", group.equal=c("loadings", "intercepts"))
cf1mg.ei <- summary(fit1mg.ei, fit.measures = TRUE, standardized = TRUE)
cf2 <- compareFit(fit1mg.el, fit1mg.ei)
t1 <- rbind(t1, cfxtract(cf2, "scalar invariance")[2,])

t1$Model[1] <- "configural invariance"
qflextable(t1)

Model	chisq	df	cfi	tli	rmsea	rmsea.ci.lo	rmsea.ci.up	srmr	d.cfi	d.tli	d.rmsea	d.srmr
configural invariance	318.01	232	0.977	0.973	0.043	0.030	0.054	0.059
metric invariance	315.33	246	0.976	0.973	0.042	0.027	0.056	0.064	-0.001	0.000	0.000	0.005
scalar invariance	339.78	260	0.972	0.971	0.045	0.030	0.057	0.067	-0.004	-0.003	0.002	0.003

Gender invariance

fit1mgs <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2:3 & sugu %in% 1:2), estimator = "WLSMV", group="sugu")
sf1mgs <- summary(fit1mg, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg)

#### Gender metric invariance 
fit1mgs.el <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2:3 & sugu %in% 1:2), estimator = "WLSMV", group="sugu", group.equal="loadings")
sf1mgs.el <- summary(fit1mgs.el, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg.el)
cf1s  <- compareFit(fit1mgs, fit1mgs.el)
t1s <- cfxtract(cf1s, "metric invariance") 


#### Gender scalar invariance
fit1mgs.ei <- cfa(mud1, data=subset(Di, vanus_grupp %in% 2:3 & sugu %in% 1:2), estimator = "WLSMV", group="sugu", group.equal=c("loadings", "intercepts"))
sf1mgs.ei <- summary(fit1mgs.ei, fit.measures = TRUE, standardized = TRUE)
cf2s <- compareFit(fit1mgs.el, fit1mgs.ei)
t1s <- rbind(t1s, cfxtract(cf2s, "scalar invariance")[2,])
qflextable(t1s)

Model	chisq	df	cfi	tli	rmsea	rmsea.ci.lo	rmsea.ci.up	srmr	d.cfi	d.tli	d.rmsea	d.srmr
	290.32	232	1.000	1.019	0.034	0.019	0.046	0.054
metric invariance	328.27	246	0.995	0.995	0.047	0.032	0.059	0.066	-0.005	-0.024	0.012	0.011
scalar invariance	373.23	260	0.986	0.985	0.053	0.041	0.065	0.071	-0.010	-0.010	0.007	0.006

Internal consistency and mean differences

Internal consistency and interitem correlations

Dis <- subset(Di, select=psc_l)
t2 <- rbind(alphaomega(Dis, att), alphaomega(Dis, int), alphaomega(Dis, ext)) |> round(2) |> as.data.frame()
t2 <- data.frame(Subscale=c("ATT", "INT", "EXT"), t2)
qflextable(t2)

Subscale	alpha	omega.tot	avgcor	mincor	maxcor
ATT	0.76	0.82	0.39	0.24	0.64
INT	0.82	0.85	0.47	0.33	0.60
EXT	0.65	0.74	0.21	0.12	0.35

Internal consistency and interitem correlations: boys

Dis <- subset(Di, sugu==1, select=psc_l)
t2 <- rbind(alphaomega(Dis, att), alphaomega(Dis, int), alphaomega(Dis, ext)) |> round(2) |> as.data.frame()
t2 <- data.frame(Subscale=c("ATT", "INT", "EXT"), t2)
qflextable(t2)

Subscale	alpha	omega.tot	avgcor	mincor	maxcor
ATT	0.78	0.83	0.41	0.25	0.61
INT	0.81	0.84	0.45	0.32	0.57
EXT	0.60	0.73	0.18	0.03	0.37

Internal consistency and interitem correlations: girls

Dis <- subset(Di, sugu==2, select=psc_l)
t2 <- rbind(alphaomega(Dis, att), alphaomega(Dis, int), alphaomega(Dis, ext)) |> round(2) |> as.data.frame()
t2 <- data.frame(Subscale=c("ATT", "INT", "EXT"), t2)
qflextable(t2)

Subscale	alpha	omega.tot	avgcor	mincor	maxcor
ATT	0.74	0.80	0.36	0.17	0.65
INT	0.82	0.86	0.48	0.36	0.60
EXT	0.69	0.77	0.25	0.08	0.42

Gender means and differences

M1 - boys M2 - girls

rbind(msdd(Di$att, Di$sugu2, "ATT"),
      msdd(Di$int, Di$sugu2, "INT"),
      msdd(Di$ext, Di$sugu2, "EXT")) |> qflextable()

Subscale	M1	SD1	M2	SD2	d	d.lo	d.up
ATT	3.97	2.45	4.72	2.47	-0.31	-0.51	-0.10
INT	3.47	2.66	4.93	2.66	-0.55	-0.76	-0.35
EXT	3.12	2.13	2.52	2.12	0.28	0.08	0.48

Age group means and differences

M1 - 11-14 M2 - 15+

rbind(msdd(Di$att, Di$vg2, "ATT"),
      msdd(Di$int, Di$vg2, "INT"),
      msdd(Di$ext, Di$vg2, "EXT")) |> qflextable()

Subscale	M1	SD1	M2	SD2	d	d.lo	d.up
ATT	3.94	2.60	4.87	2.29	-0.38	-0.58	-0.17
INT	3.55	2.66	4.97	2.68	-0.53	-0.74	-0.33
EXT	2.72	2.12	2.87	2.16	-0.07	-0.27	0.13

Correlations with RCADS and Kidscreen

vars <- c("att", "int", "ext", "RCADS.Depressioon", "RCADS.Ärevus", "RCADS.kokku", "Kidscreen.HEAOLU")
kt <- psych::corr.test(subset(Di, vg2%in%2:3)[, vars])
ktr <- kt$r
ktr[upper.tri(ktr)]<-NA
rownames(ktr) <- colnames(ktr) <- c("att", "int", "ext", "RCADS.Depr", "RCADS.Anx", "RCADS.tot", "KS.WB")
ktr <- data.frame(Scale = rownames(ktr), round(ktr,3))
maxp <- format(max(kt$p), digits=3, scientific=FALSE)
qflextable(ktr)

Scale	att	int	ext	RCADS.Depr	RCADS.Anx	RCADS.tot	KS.WB
att	1.000
int	0.615	1.000
ext	0.478	0.342	1.000
RCADS.Depr	0.679	0.785	0.429	1.000
RCADS.Anx	0.638	0.747	0.354	0.788	1.000
RCADS.tot	0.693	0.808	0.410	0.933	0.957	1.000
KS.WB	-0.449	-0.592	-0.237	-0.647	-0.513	-0.605	1

The largest p-value for these correlations was 0.00000229. N varies from 386 to 390.

Parent reports :CFA

mud2 <- " ATT =~ LV6 + LV7 + LV8 + LV9 + LV10
         INT =~ LV1 + LV2 + LV3 + LV4 + LV5
         EXT =~ LV11 + LV12 + LV13 + LV14 + LV15 + LV16 + LV17
"

fit2 <- cfa(mud2, data=subset(Dj, vg2 %in% 2:3), estimator="WLSMV")

Warning in lav_options_set(opt): lavaan WARNING: estimator "WLSMV" is not recommended for continuous data.
          Did you forget to set the ordered= argument?

getfms(fit2)

chisq.sc	df.sc	pvalue.sc	cfi.r	tli.r	srmr	rmsea.r	rmsea.ci.lo.r	rmsea.ci.up.r
194	116	0.0000082	0.995	0.995	0.0647	0.0475	0.0355	0.0591

semPaths(fit2, rotation=2, what="std", residuals=FALSE, edge.color="black", esize=1, edge.label.cex=1, sizeMan = 6, node.width=1, node.height=1, weighted=FALSE)

ss <- standardizedsolution(fit2) |> as.data.frame()
ss[,] <- lapply(ss, formatif, 3)
qflextable(ss)

lhs	op	rhs	est.std	se	z	pvalue	ci.lower	ci.upper
ATT	=~	LV6	0.764	0.0456	16.77	0.000000000000000000	0.675	0.854
ATT	=~	LV7	0.325	0.0767	4.23	0.000023120224550865	0.174	0.475
ATT	=~	LV8	0.757	0.0361	20.97	0.000000000000000000	0.686	0.828
ATT	=~	LV9	0.786	0.0436	18.05	0.000000000000000000	0.701	0.872
ATT	=~	LV10	0.587	0.0598	9.81	0.000000000000000000	0.470	0.704
INT	=~	LV1	0.700	0.0589	11.88	0.000000000000000000	0.585	0.816
INT	=~	LV2	0.820	0.0460	17.82	0.000000000000000000	0.730	0.910
INT	=~	LV3	0.553	0.0539	10.26	0.000000000000000000	0.447	0.659
INT	=~	LV4	0.471	0.0709	6.64	0.000000000032114311	0.332	0.609
INT	=~	LV5	0.669	0.0581	11.52	0.000000000000000000	0.555	0.783
EXT	=~	LV11	0.595	0.0660	9.02	0.000000000000000000	0.466	0.724
EXT	=~	LV12	0.705	0.0498	14.15	0.000000000000000000	0.607	0.802
EXT	=~	LV13	0.633	0.0527	12.01	0.000000000000000000	0.530	0.736
EXT	=~	LV14	0.702	0.0512	13.71	0.000000000000000000	0.602	0.803
EXT	=~	LV15	0.714	0.0399	17.90	0.000000000000000000	0.636	0.792
EXT	=~	LV16	0.604	0.0542	11.15	0.000000000000000000	0.498	0.710
EXT	=~	LV17	0.546	0.0677	8.06	0.000000000000000666	0.413	0.678
LV6	~~	LV6	0.416	0.0696	5.97	0.000000002337326332	0.279	0.552
LV7	~~	LV7	0.895	0.0498	17.96	0.000000000000000000	0.797	0.992
LV8	~~	LV8	0.427	0.0547	7.80	0.000000000000005995	0.319	0.534
LV9	~~	LV9	0.381	0.0685	5.57	0.000000026144011844	0.247	0.516
LV10	~~	LV10	0.655	0.0703	9.33	0.000000000000000000	0.518	0.793
LV1	~~	LV1	0.510	0.0825	6.18	0.000000000633316288	0.348	0.672
LV2	~~	LV2	0.327	0.0755	4.33	0.000014704537161903	0.179	0.475
LV3	~~	LV3	0.694	0.0596	11.65	0.000000000000000000	0.577	0.811
LV4	~~	LV4	0.779	0.0667	11.67	0.000000000000000000	0.648	0.909
LV5	~~	LV5	0.553	0.0776	7.12	0.000000000001084022	0.401	0.705
LV11	~~	LV11	0.646	0.0785	8.23	0.000000000000000222	0.492	0.800
LV12	~~	LV12	0.504	0.0702	7.18	0.000000000000711653	0.366	0.641
LV13	~~	LV13	0.599	0.0667	8.98	0.000000000000000000	0.468	0.730
LV14	~~	LV14	0.506	0.0720	7.04	0.000000000001975753	0.365	0.648
LV15	~~	LV15	0.490	0.0569	8.61	0.000000000000000000	0.379	0.602
LV16	~~	LV16	0.635	0.0654	9.72	0.000000000000000000	0.507	0.764
LV17	~~	LV17	0.702	0.0739	9.50	0.000000000000000000	0.557	0.847
ATT	~~	ATT	1.000	0.0000	NA	NA	1.000	1.000
INT	~~	INT	1.000	0.0000	NA	NA	1.000	1.000
EXT	~~	EXT	1.000	0.0000	NA	NA	1.000	1.000
ATT	~~	INT	0.362	0.0812	4.46	0.000008173817187318	0.203	0.521
ATT	~~	EXT	0.543	0.0591	9.20	0.000000000000000000	0.427	0.659
INT	~~	EXT	0.299	0.0804	3.71	0.000205557733914397	0.141	0.456

Age group invariance

fit2mg <- cfa(mud2, data=subset(Dj, vanus_grupp %in% 1:3), estimator = "WLSMV", group="vanus_grupp")
sf2mg <- summary(fit2mg, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg)

#### Age group metric invariance (multigroup model, equal loadings)
fit2mg.el <- cfa(mud2, data=subset(Dj, vanus_grupp %in% 1:3), estimator = "WLSMV", group="vanus_grupp", group.equal="loadings")
sf2mg.el <- summary(fit2mg.el, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg.el)
cf3  <- compareFit(fit2mg, fit2mg.el)
t3 <- cfxtract(cf3, "metric invariance") 
anova(fit2mg, fit2mg.el)

Scaled Chi-Squared Difference Test (method = "satorra.2000")

lavaan NOTE:
    The "Chisq" column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
           Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)
fit2mg    348         291.98                              
fit2mg.el 376         356.85     28.087      28     0.4599

#### Age group scalar invariance (multigroup model, equal intercepts) 
fit2mg.ei <- cfa(mud2, data=subset(Dj, vanus_grupp %in% 1:3), estimator = "WLSMV", group="vanus_grupp", group.equal=c("loadings", "intercepts"))
cf2mg.ei <- summary(fit2mg.ei, fit.measures = TRUE, standardized = TRUE)
cf4 <- compareFit(fit2mg.el, fit2mg.ei)
anova(fit2mg.el, fit2mg.ei)

Scaled Chi-Squared Difference Test (method = "satorra.2000")

lavaan NOTE:
    The "Chisq" column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
           Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)
fit2mg.el 376         356.85                              
fit2mg.ei 404         376.79     31.959      28     0.2762

t4 <- rbind(t3, cfxtract(cf4, "scalar invariance")[2,])

t4$Model[1] <- "configural invariance"
qflextable(t4)

Model	chisq	df	cfi	tli	rmsea	rmsea.ci.lo	rmsea.ci.up	srmr	d.cfi	d.tli	d.rmsea	d.srmr
configural invariance	465.12	348	0.963	0.956	0.053	0.039	0.065	0.075
metric invariance	458.10	376	0.963	0.960	0.051	0.032	0.066	0.082	0.000	0.003	-0.002	0.007
scalar invariance	488.09	404	0.962	0.962	0.049	0.031	0.065	0.084	-0.001	0.002	-0.001	0.002

Gender invariance

fit2mgs <- cfa(mud2, data=subset(Dj, vanus_grupp %in% 1:3 & sugu %in% 1:2), estimator = "WLSMV", group="sugu")
sf2mgs <- summary(fit2mg, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg)

#### Gender metric invariance 
fit2mgs.el <- cfa(mud2, data=subset(Dj, vanus_grupp %in% 1:3 & sugu %in% 1:2), estimator = "WLSMV", group="sugu", group.equal="loadings")
sf2mgs.el <- summary(fit2mgs.el, fit.measures = TRUE, standardized = TRUE)
#getfms(sf1mg.el)
cf2s  <- compareFit(fit2mgs, fit2mgs.el)
anova(fit2mgs, fit2mgs.el)

Scaled Chi-Squared Difference Test (method = "satorra.2000")

lavaan NOTE:
    The "Chisq" column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
            Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)
fit2mgs    232         204.90                              
fit2mgs.el 246         240.01      14.94      14     0.3823

t2s <- cfxtract(cf2s, "metric invariance") 


#### Gender scalar invariance
fit2mgs.ei <- cfa(mud2, data=subset(Dj, vanus_grupp %in% 1:3 & sugu %in% 1:2), estimator = "WLSMV", group="sugu", group.equal=c("loadings", "intercepts"))
sf2mgs.ei <- summary(fit2mgs.ei, fit.measures = TRUE, standardized = TRUE)
cf2s <- compareFit(fit2mgs.el, fit2mgs.ei)
anova (fit2mgs.el, fit2mgs.ei)

Scaled Chi-Squared Difference Test (method = "satorra.2000")

lavaan NOTE:
    The "Chisq" column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
            Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
fit2mgs.el 246         240.01                                  
fit2mgs.ei 260         274.25     52.414      14  2.392e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

t2s <- rbind(t2s, cfxtract(cf2s, "scalar invariance")[2,])
qflextable(t2s)

Model	chisq	df	cfi	tli	rmsea	rmsea.ci.lo	rmsea.ci.up	srmr	d.cfi	d.tli	d.rmsea	d.srmr
	339.27	232	0.971	0.966	0.047	0.036	0.057	0.064
metric invariance	327.32	246	0.970	0.967	0.046	0.032	0.059	0.069	-0.001	0.001	-0.001	0.005
scalar invariance	363.01	260	0.962	0.960	0.051	0.038	0.063	0.073	-0.008	-0.007	0.004	0.004

Parental reports: Correlations with RCADS and Kidscreen

Age groups 2 and 3

vars <- c("attlv", "intlv", "extlv", "RCADS.Depressioon", "RCADS.Ärevus", "RCADS.kokku", "Kidscreen.HEAOLU")
kt <- psych::corr.test(subset(Dij, vg2%in%2:3)[, vars])
ktr <- kt$r
ktr[upper.tri(ktr)]<-NA
rownames(ktr) <- colnames(ktr) <- c("att", "int", "ext", "RCADS.Depr", "RCADS.Anx", "RCADS.tot", "KS.WB")
p <- kt$p[, 1:3]
p[upper.tri(p)] <- NA
ktr <- data.frame(Scale = rownames(ktr), do.call(cbind, lapply(1:3, \(i) data.frame(round(ktr[,i],2), round(p[,i], 5)))))
names(ktr)[-1] <- c("att.r", "att.p", "int.r", "int.p", "ext.r", "ext.p")
qflextable(ktr)

Scale	att.r	att.p	int.r	int.p	ext.r	ext.p
att	1.00	0.00000
int	0.25	0.00007	1.00	0.00000
ext	0.44	0.00000	0.21	0.00061	1.00	0.00000
RCADS.Depr	0.23	0.00017	0.34	0.00000	0.12	0.06651
RCADS.Anx	0.16	0.00865	0.31	0.00000	0.09	0.14727
RCADS.tot	0.21	0.00094	0.34	0.00000	0.11	0.08707
KS.WB	-0.15	0.01959	-0.32	0.00000	-0.13	0.04018

n <- kt$n[,1:3]

N varies from 254 to 255.

Age group 1

vars <- c("attlv", "intlv", "extlv", "Kidscreen.HEAOLU", "Kidscreen.MEELEOLU")
kt <- psych::corr.test(subset(Dij, vanus_grupp%in%1)[, vars])
ktr <- kt$r
ktr[upper.tri(ktr)]<-NA
rownames(ktr) <- colnames(ktr) <- c("att", "int", "ext", "KS.WB", "KS.Mood")
p <- kt$p[, 1:3]
p[upper.tri(p)] <- NA
ktr <- data.frame(Scale = rownames(ktr), do.call(cbind, lapply(1:3, \(i) data.frame(round(ktr[,i],2), round(p[,i], 5)))))
names(ktr)[-1] <- c("att.r", "att.p", "int.r", "int.p", "ext.r", "ext.p")
qflextable(ktr)

Scale	att.r	att.p	int.r	int.p	ext.r	ext.p
att	1.00	0.00000
int	0.40	0.00001	1.00	0.00000
ext	0.57	0.00000	0.41	0.00001	1.00	0.00000
KS.WB	-0.14	0.13086	-0.21	0.03099	-0.14	0.14727
KS.Mood	-0.22	0.02205	-0.15	0.11896	-0.17	0.08261

n <- kt$n[,1:3]

N varies from 108 to 111.

Correlations between self- and parental reports

Included: only respondents in 2 older age groups if both child and parent responded in Estonian

#andmetabeli mõõtmed, st read x veerud
dim(subset(Dij, vg2%in% 2:3))

[1] 384 473

ct <- corr.test(Dij[, c("att", "int", "ext")], Dij[, c("attlv", "intlv", "extlv")])
res <- data.frame(var = rownames(ct$r), r = round(diag(ct$r),3), p = round(diag(ct$p), 5), N = if(length(ct$n)==1) ct$n else diag(ct$n))
qflextable(res)

var	r	p	N
att	0.272	0.00001	254
int	0.344	0.00000	254
ext	0.271	0.00001	254

Correlations between single items (the order of items was different in child’s and parental questionnaires: in the table, the order used in child’s version is used).

pk <- openxlsx::read.xlsx("../LVTU_tunnuste_kirjeldused.xlsx", sheet="PSC")
ct <- corr.test(Dij[, pk$laps], Dij[, pk$lv])
res <- data.frame(var = rownames(ct$r), item = pk$laps_item, r = round(diag(ct$r),3), p = round(diag(ct$p), 5), N = if(length(ct$n)==1) ct$n else diag(ct$n))
qflextable(res)

var	item	r	p	N
PSC1	1. Olen rahutu, ei suuda paigal püsida	0.324	0.00000	255
PSC2	2. Olen kurb, õnnetu	0.168	0.00717	255
PSC3	3. Unistan liiga palju	0.073	0.24758	254
PSC4	4. Keeldun oma asju jagamast	0.170	0.00670	254
PSC5	5. Ma ei mõista teiste inimeste tundeid	0.141	0.02442	255
PSC6	6. Tunnen lootusetust	0.327	0.00000	250
PSC7	7. Mul on raske keskenduda	0.324	0.00000	253
PSC8	8. Kaklen teiste lastega	0.117	0.06290	254
PSC9	9. Olen enda suhtes kriitiline, karm	0.252	0.00005	254
PSC10	10. Süüdistan oma muredes teisi	0.233	0.00018	253
PSC11	11. Mul ei ole nii lõbus kui varem	0.156	0.01290	254
PSC12	12. Ma ei järgi reegleid	0.218	0.00049	252
PSC13	13. Tunnen vajadust pidevalt tegutseda, nagu oleksin üles keeratud	0.105	0.09414	255
PSC14	14. Narrin teisi	0.186	0.00306	253
PSC15	15. Muretsen palju	0.234	0.00017	254
PSC16	16. Võtan asju, mis ei kuulu mulle	0.223	0.00034	254
PSC17	17. Minu tähelepanu kaldub kergesti kõrvale	0.277	0.00001	255

Prediction of health problems

Parental questionnaire included a question about different diseases / disorders diagnosed in the child; among these: anxiety disorders, mood disorders, eating disorders, learning difficulties, ADHD, autism spectrum disorders.

• LV_terv_prob_5 343 13.5. …Ärevushäire, nt mõni foobia, obsessiiv-kompulsiivne häire või paanikahäire
• LV_terv_prob_6 344 13.6. …Meeleoluhäire, nt depressioon, bipolaarne häire, maania või düstüümia
• LV_terv_prob_7 345 13.7. …Söömishäire, nt anoreksia, buliimia või liigsöömishäire
• LV_terv_prob_8 346 13.8. …Õpiraskused või spetsiifilise akadeemilise oskuse kasutamise raskused (nt kõneareng, lugemine, kirjutamine, aritmeetika)
• LV_terv_prob_9 347 13.9. …Aktiivsus- ja tähelepanuhäire ehk ATH
• LV_terv_prob_10 348 13.10. …Autismispektri häire (varasemalt nimetatud ka kui autism, Aspergeri sündroom, pervasiivne arenguhäire)

Descriptive stats about health problems

Table: y1= cases of this diagnosis in the first age group; n1 = cases without this diagnosis in the first age group … etc

vars <- paste0("LV_terv_prob_", 5:10)
legend <- c("anxiety disorders", "mood disorders", "eating disorders", "learning difficulties", "ADHD", "autism spectrum disorders")
res <- data.frame(t(aggregate(Dij[, vars], Dij['vanus_grupp'], sum, na.rm=T)[,-1]))
res2 <- data.frame(t(aggregate(Dij[, vars], Dij['vanus_grupp'], function(x) sum(1-x, na.rm=T))[,-1]))
RES <- do.call(data.frame,mapply(data.frame, res, res2 , SIMPLIFY = FALSE))
names(RES) <- c(outer(c("y", "n"), 1:3, FUN=paste0))
RES$Y = rowSums(res)
RES$N = rowSums(res2)

RES <- data.frame(diagn=legend, RES)
qflextable(RES)

diagn	y1	n1	y2	n2	y3	n3	Y	N
anxiety disorders	5	106	3	131	11	111	19	348
mood disorders	2	109	3	131	17	105	22	345
eating disorders	4	107	2	132	2	120	8	359
learning difficulties	11	100	11	123	6	116	28	339
ADHD	5	106	8	126	6	116	19	348
autism spectrum disorders	2	109	4	130	4	118	10	357

Predicting internalizing disorders (anxiety + mood disorders)

Just sex and age as predictors (using only age groups 2 and 3!)

glm.anx1 <- glm(I(LV_terv_prob_5 | LV_terv_prob_6)~ factor(sugu2) + vanus, data=subset(Dij, vg2 %in% 2:3), family="binomial")
summary(glm.anx1)

Call:
glm(formula = I(LV_terv_prob_5 | LV_terv_prob_6) ~ factor(sugu2) + 
    vanus, family = "binomial", data = subset(Dij, vg2 %in% 2:3))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -7.8256     1.7534  -4.463 8.08e-06 ***
factor(sugu2)2   0.6950     0.4556   1.525  0.12715    
vanus            0.3532     0.1115   3.167  0.00154 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 167.11  on 250  degrees of freedom
Residual deviance: 151.99  on 248  degrees of freedom
  (133 observations deleted due to missingness)
AIC: 157.99

Number of Fisher Scoring iterations: 5

PseudoR2(glm.anx1)

 McFadden 
0.0904597 

Adding child’s self-reported PSC-internalizing scale

glm.anx2 <- glm(I(LV_terv_prob_5 | LV_terv_prob_6)~  factor(sugu2) + vanus + scale(int), data=subset(Dij, vg2 %in% 2:3), family="binomial") 
summary(glm.anx2)

Call:
glm(formula = I(LV_terv_prob_5 | LV_terv_prob_6) ~ factor(sugu2) + 
    vanus + scale(int), family = "binomial", data = subset(Dij, 
    vg2 %in% 2:3))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -7.0928     1.7861  -3.971 7.15e-05 ***
factor(sugu2)2   0.6076     0.4617   1.316   0.1882    
vanus            0.3002     0.1142   2.628   0.0086 ** 
scale(int)       0.4277     0.2267   1.887   0.0592 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 166.89  on 249  degrees of freedom
Residual deviance: 148.18  on 246  degrees of freedom
  (134 observations deleted due to missingness)
AIC: 156.18

Number of Fisher Scoring iterations: 6

PseudoR2(glm.anx2)

 McFadden 
0.1121497 

–> OR = 1.67 (for 1 SD change in child’s self-report)

Adding parent report

glm.anx3 <- glm(I(LV_terv_prob_5 | LV_terv_prob_6)~  factor(sugu2) + vanus + scale(int) + scale(intlv), data=subset(Dij, vg2 %in% 2:3), family="binomial") 
summary(glm.anx3)

Call:
glm(formula = I(LV_terv_prob_5 | LV_terv_prob_6) ~ factor(sugu2) + 
    vanus + scale(int) + scale(intlv), family = "binomial", data = subset(Dij, 
    vg2 %in% 2:3))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -8.5034     2.0265  -4.196 2.72e-05 ***
factor(sugu2)2   0.2315     0.4978   0.465 0.641883    
vanus            0.3788     0.1266   2.993 0.002762 ** 
scale(int)       0.2463     0.2430   1.014 0.310807    
scale(intlv)     1.0285     0.2699   3.811 0.000139 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 162.33  on 248  degrees of freedom
Residual deviance: 124.41  on 244  degrees of freedom
  (135 observations deleted due to missingness)
AIC: 134.41

Number of Fisher Scoring iterations: 6

PseudoR2(glm.anx3)

 McFadden 
0.2335861 

–> OR = 3.15 (for 1 SD change in parent report)

Predicting ADHD from attention subscale

glm.ath1 <- glm(LV_terv_prob_9 ~ factor(sugu2) + vanus, data=subset(Dij, vg2 %in% 2:3), family="binomial") 
summary(glm.ath1)

Call:
glm(formula = LV_terv_prob_9 ~ factor(sugu2) + vanus, family = "binomial", 
    data = subset(Dij, vg2 %in% 2:3))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)
(Intercept)     -1.0086     1.7940  -0.562    0.574
factor(sugu2)2  -0.3457     0.5590  -0.618    0.536
vanus           -0.1188     0.1301  -0.913    0.361

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 108.02  on 250  degrees of freedom
Residual deviance: 106.65  on 248  degrees of freedom
  (133 observations deleted due to missingness)
AIC: 112.65

Number of Fisher Scoring iterations: 6

PseudoR2(glm.ath1)

  McFadden 
0.01269651 

Adding self-report

glm.ath2 <- glm(LV_terv_prob_9 ~ factor(sugu2) + vanus + scale(att), data=subset(Dij, vg2 %in% 2:3), family="binomial") 
summary(glm.ath2)

Call:
glm(formula = LV_terv_prob_9 ~ factor(sugu2) + vanus + scale(att), 
    family = "binomial", data = subset(Dij, vg2 %in% 2:3))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)  
(Intercept)     -0.4205     1.8448  -0.228    0.820  
factor(sugu2)2  -0.3967     0.5650  -0.702    0.483  
vanus           -0.1641     0.1349  -1.216    0.224  
scale(att)       0.4871     0.2708   1.799    0.072 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.91  on 249  degrees of freedom
Residual deviance: 103.26  on 246  degrees of freedom
  (134 observations deleted due to missingness)
AIC: 111.26

Number of Fisher Scoring iterations: 6

PseudoR2(glm.ath2)

  McFadden 
0.04304296 

Adding parent report

glm.ath3 <- glm(LV_terv_prob_9 ~ factor(sugu2) + vanus + scale(att) + scale(attlv), data=subset(Dij, vg2 %in% 2:3), family="binomial") 
summary(glm.ath3)

Call:
glm(formula = LV_terv_prob_9 ~ factor(sugu2) + vanus + scale(att) + 
    scale(attlv), family = "binomial", data = subset(Dij, vg2 %in% 
    2:3))

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -3.86399    2.30009  -1.680    0.093 .  
factor(sugu2)2 -0.53616    0.63236  -0.848    0.397    
vanus           0.01926    0.15530   0.124    0.901    
scale(att)      0.07749    0.30639   0.253    0.800    
scale(attlv)    1.48335    0.33447   4.435 9.21e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.793  on 248  degrees of freedom
Residual deviance:  77.497  on 244  degrees of freedom
  (135 observations deleted due to missingness)
AIC: 87.497

Number of Fisher Scoring iterations: 7

PseudoR2(glm.ath3)

 McFadden 
0.2810518