RVTU ülesöömise väidete IRT analüüs 29.05.2025

Kokkuvõte: meil on andmed koduse keele kohta ainult 1. lainest. Vastamiskeel on usaldusväärselt teada ainult pabertesti vastajate kohta (“postivalim”), kuna limesurvey ei registreerinud seda täpselt (vastaja võis testi alguses vastamiskeelt vahetada ning seda ei ole registreeritud). Seega määrame tõenäolise vastamiskeele kaudselt:

  1. Kõigepealt võtame aluseks rahvastikuregistris märgitud emakeele
  2. Vastajatel, kellel on rahvastikuregistris emakeel puudu, võtame aluseks RR rahvuse (eestlastele -> eesti keel, venelastele -> vene keel, teistel jääb tühjaks)
  3. Vastajatel, kellel on 1. laines märgitud “peamine kodune keel”, kasutame seda RR emakeele asemel
  4. Vastajad, kellel on 1. laines märgitud “muu” kodune keel, märgime vastamiskeeleks eesti, kui tekstivastusest selgub, et eesti keel on tema esimene keel ning ta on vastanud eesti keeles (nt “võro”, või “eesti ja vene” või “inglise, elukaaslasega suhtlemiseks”). Kui koduseks keeleks on märgitud eesti ja vene, aga vastus on kirjutatud vene keeles, siis märgime vastamiskeeleks vene keele.

Selliselt muudetud vastamiskeele tunnus on lisatud andmetabelile (tunnus L$keel). Järgnevates analüüsides on kasutatud ainult ridu, mille puhul vastamiskeele väärtus on EESTI (potentsiaalselt 15540, kuid neist üle poole ei ole vastanud küsitlusele).

Kood 
table(L$keel)

EESTI  VENE 
15540  2662

Ülesöömise väited

Kood 
vaited <- key(L, NAME %in% paste0("L1EAT",1:5))[, c("NAME", "LAB", "LEVELS")]
vaited$LAB <- gsub("23. Järgnevalt soovime teada saada, kui suurt rolli on mänginud toit ja söömine Teie elus viimase kolme (3) kuu jooksul.","", vaited$LAB, fixed=TRUE)
qflextable(vaited)

NAME

LAB

LEVELS

L1EAT1

[Kas Te olete veetnud märkimisväärselt palju aega kehakaalust ja toidust mõeldes]

1:ei|2:jah

L1EAT2

[Kas Te olete viimase kolme kuu jooksul oma toitumist märkimisväärselt piiranud]

1:ei|2:jah

L1EAT3

[Kas Teil on viimase kolme kuu jooksul esinenud söömishoogusid, mil sõite tavalisest suurema toidukoguse]

1:ei|2:jah

L1EAT4

Kas nende söömishoogude ajal olete tundnud, et Te ei suuda oma söömist kontrollida?                                                                                      

1:ei|2:jah

L1EAT5

Kas olete viimase kolme kuu jooksul oma kehakaalu kontrollimiseks tahtlikult oksendanud, kasutanud lahtisteid või söögiisu vähendavaid preparaate?

1:ei|2:jah

Eeltöötlus

Kood 
eat1 <- subset(L, incl)[, paste0("L1EAT", 1:5)]
names(eat1) <- paste0("e", 1:5)
keep1 <- rowSums(!is.na(eat1))>0
eat1x <- eat1[keep1,]

eat1$e34 <- ifelse(eat1$e3 %in% 1, eat1$e3, eat1$e4+1)
eat1$e4[!is.na(eat1$e3) & !eat1$e4 %in% 1:2] <- 1
eat1a <- eat1[, 1:5]
eat1 <- eat1[, c("e1", "e2", "e34", "e5")]


eat2 <- subset(L, incl)[, paste0("L2EAT", 1:5)]
names(eat2) <- paste0("e", 1:5)
eat2$e34 <- ifelse(eat2$e3 %in% 1, eat2$e3, eat2$e4+1)
eat2$e4[!is.na(eat2$e3) & !eat2$e4 %in% 1:2] <- 1
eat2a <- eat2[, 1:5]
eat2 <- eat2[, c("e1", "e2", "e34", "e5")]

eat3 <- subset(L, incl)[, paste0("L3EAT", 1:5)]
names(eat3) <- paste0("e", 1:5)
eat3$e34 <- ifelse(eat3$e3 %in% 1, eat3$e3, eat3$e4+1)
eat3$e4[!is.na(eat3$e3) & !eat3$e4 %in% 1:2] <- 1
eat3a <- eat3[, 1:5]
eat3 <- eat3[, c("e1", "e2", "e34", "e5")]

Enne analüüse tuleb kuidagi toime tulla 3. ja 4. väite sõltuvusega: 4. väide esitati ainult neile, kes vastasid 3. väitele jaatavalt. Lähtume eeldusest, et kõik, kes vastasid 3. küsimusele eitavalt, oleksid vastanud eitavalt ka 4. küsimusele - st vastajatel, kellel ei ole esinenud söömishooge, ei ole esinenud ka kontrollikaotusega söömishooge. Seega asendame kõigil, kes 3. küsimusele vastasid eitavalt, 4. küsimuse väärtuse NA väärtusega 1 (= ei).

Kood 
table(uus=eat1a$e4, vana=subset(L, incl)$L1EAT4, exclude=NULL)[1:2,]  |> unclass() |> as.data.frame() |> flexit()

X.

X1

X2

NA.

1

591

0

3 625

2

0

647

0

Teine võimalus on kombineerida 3. ja 4. väide üheks tunnuseks, kus

  • 1 = ei ole esinenud söömishoogusid
  • 2 = on esinenud söömishooge, kuid ilma kontrollikaotuseta
  • 3 = on esinenud söömishooge ja kontrolli kaotamise tunne
Kood 
t1 <- unclass(table(e34=eat1$e34, e3=subset(L, incl)$L1EAT3))
t2 <- unclass(table(e34=eat1$e34, e4=subset(L, incl)$L1EAT4, exclude=NULL)[1:3,])
tbl <- data.frame(E34=1:3, cbind(t1, NA, t2))
names(tbl) <- c("e34", "e3=1", "e3=2", " ", "e4=1", "e4=2", "e4<NA>")
qflextable(tbl)

e34

e3=1

e3=2

e4=1

e4=2

e4<NA>

1

3 621

0

0

0

3 621

2

0

591

591

0

0

3

0

647

0

647

0

Eeldused jm 1. laine andmetel

Dimensionaalsus ja IRT mudeli valik

Võrdleme dimensionaalsuse kontrolliks 1- ja 2-mõõtmelist mirt (2PL) mudelit:

Kood 
load("2025-05-29-precomputed.rda")
eat1ax <- na.omit(eat1)
# m1 <- mirt(eat1ax, model=1, itemtype="graded", calcNull=TRUE)
# m2 <- mirt(eat1ax, model=2, itemtype="graded", calcNull=TRUE, technical  =list(NCYCLES=c(10000, 2000)))
# Iteration: 2559, Log-Lik: -9582.269, Max-Change: 0.00008
anova(m1, m2) |> flexit()

X.

AIC

SABIC

HQ

BIC

logLik

X2

df

p

m1

16 996.44

17 026.19

17 016.92

17 054.79

-8 489.218

m2

16 948.76

16 988.43

16 976.08

17 026.56

-8 462.380

53.675

3

0

Andmed on selgelt vähemalt kahemõõtmelised (kahefaktoriline eksploratiivne mudel sobitub andmetega paremini kui ühemõõtmeline).

Eksploratiivse kahemõõtmelise mudeli puhul on faktorstruktuur selline:

Kood 
smm2 <- summary(m2, verbose=FALSE)
cbind(smm2$rotF, smm2$h2) |> flexit()

X.

F1

F2

h2

e1

0.575

0.379

0.688

e2

0.955

-0.095

0.833

e34

-0.005

0.959

0.915

e5

0.556

0.301

0.564

Korrelatsioon faktorite vahel (oblimin) r= 0.49.

Mudelite sobitusastme näitajad:

Kood 
data.frame(m=c("One factor", "Two factors"), rbind(fitses(m1,3), fitses(m2,3))) |> qflextable()

m

RMSEA

CFI

TLI

BIC

AIC

One factor

0.036

0.950

0.936

17 054.79

16 996.44

Two factors

0.026

0.979

0.966

17 026.56

16 948.76

Kuigi kahefaktoriline mudel on “oluliselt” parem, siis on ka ühefaktorilise mudeli sobitusastme näitajad täiesti aktsepteeritavad. Lisaks: eesmärk on eelkõige leida 5 söömispatoloogia küsimuse ühisosa, mitte kõigi nüansside modelleerimine.

Sisuliselt samasuguse mudeli saab teha andmestikuga, kus on säilitatud kõik 5 küsimust eraldi dihhotoomisete tunnustega ja 4. küsimuse puuduvad andmed on asendatud vastavalt eespool kirjeldatud skeemile. Sel juhul tuleb arvestada, et 3. ja 4. küsimuse vahel on meetodist (sh küsimuse esitamise viisist) tulenev sõltuvus, mida mudel ei tohiks ignoreerida. Kinnitavas faktoranalüüsis võiks seda kujutada 3. ja 4. tunnuse korreleeritud jääkidena. IRT mudelis (kasutades mirt paketti) on võimalus kujutada sisuliselt sama seost “nuditud” bifaktori mudeli abil, kus lisaks üldfaktorile (mis on seotud kõgi 5 küsimusega) on ainult üks spetsiifiline faktor, millel on laadungid ainult 3. ja 4. tunnusel. Selline bifaktori mudel on sisuliselt samaväärne mudeliga, kus on lubatud korrelatsioon 3. ja 4. tunnuse jääkide vahel (nt saranselt jääkkorrelatsiooniga on spetsiifiline faktor sõltumatu üldfaktorist).

Kood 
ldgs(mbf,3) |> qflextable()

G

F1

h2

e1

0.991

0.000

0.982

e2

0.671

0.000

0.451

e3

0.594

0.800

0.993

e4

0.685

0.726

0.996

e5

0.750

0.000

0.563

Võrdleme bifaktori mudeli puhul erinevaid IRT parametriseeringuid: ühe- ja kaheparameetrilist logistilist mudelit, alumise või ülemise asümptoodiga kolmeparameetrilist mudelit ning neljaparameetrilist mudelit.

Kood 
anova(mbf_rasch, mbf, mbf_3pl, mbf_3plu, mbf_4pl) |> flexit()

X.

AIC

SABIC

HQ

BIC

logLik

X2

df

p

mbf_rasch

17 499.68

17 522.83

17 515.62

17 545.07

-8 742.840

mbf

16 987.77

17 027.45

17 015.09

17 065.58

-8 481.884

521.912

5

0.000

mbf_3pl

16 983.59

17 039.81

17 022.29

17 093.82

-8 474.794

14.180

5

0.015

mbf_3plu

16 987.12

17 043.33

17 025.81

17 097.35

-8 476.558

-3.527

0

mbf_4pl

16 985.42

17 058.17

17 035.50

17 128.08

-8 470.709

11.698

5

0.039

Üheparameetriline mudel on selgelt ebapiisav - st ei saa eeldada, et küsimused erinevad ainult raskusastme poolest. 3- ja 4-parameetrilised mudelid on log tõepära alusel 2PL mudelist paremad, kuid BIC ja SABIC (mis “karistavad” mudelit lisaparameetrite eest rohkem kui nt AIC) on kõige väiksemad 2PL mudelil.

Lõpuks: bifaktori mudelite sobitusastme näitajad:

Kood 
data.frame(model=c("1PL", "2PL", "3PL", "3PLu", "4PL"), rbind(fitses(mbf_rasch,3), fitses(mbf,3), fitses(mbf_3pl, 3), fitses(mbf_3plu,3), fitses(mbf_4pl, 3))) |> qflextable()

model

RMSEA

CFI

TLI

BIC

AIC

1PL

0.069

0.855

0.843

17 545.07

17 499.68

2PL

0.021

0.990

0.986

17 065.58

16 987.77

3PL

0.021

0.992

0.985

17 093.82

16 983.59

3PLu

0.023

0.991

0.983

17 097.35

16 987.12

4PL

0.025

0.993

0.979

17 128.08

16 985.42

Lokaalse sõltumatuse eeldus

Chen ja Thissen (1997) G2-statistik, mis nullhüpoteesi korral (st kui lokaalset sõltuvust ei esine = lokaalse sõltumatuse eeldus on täidetud) käitub nagu χ2-jaotus 1 vabadusastmega. Seega peaksime muretsema, kui mõne väitepaari vahel on G2 suurem kui 3.8414588 (5% olulisuse nivoo korral) või 6.6348966 (kasutades 1% olulisuse nivood).

Kood 
residuals(m1, type="LDG2", verbose=FALSE, approx.z=!TRUE) |> (\(x) ("[<-"(x, upper.tri(x), NA)))() |> flexit()

X.

e1

e2

e34

e5

e1

e2

0.584

e34

36.996

15.677

e5

3.460

3.318

5.979

Kood 
residuals(mbf, type="LDG2", verbose=FALSE, approx.z=!TRUE) |> (\(x) ("[<-"(x, upper.tri(x), NA)))() |> flexit()

X.

e1

e2

e3

e4

e5

e1

e2

0.453

e3

0.586

7.988

e4

0.689

9.795

0.427

e5

3.653

3.149

0.929

1.945

See tähendab, et lokaalse sõltumatuse eeldus ei ole päriselt täidetud. Bifaktori mudeli puhul on kõige suuremad jääkkorrelatsioonid 2. ning 3. ja 4. küsimuse vahel. Seda olukorda saaks parandada nii, et lisada bifaktori mudelis ka 2. küsimus spetsiifilise faktori indikaatoriks. Sellise mudeli faktorstruktuur:

Kood 
ldgs(mbf2,3) |> flextable()

G

F1

h2

e1

0.859

0.000

0.738

e2

0.780

-0.320

0.711

e3

0.711

0.695

0.988

e4

0.766

0.641

0.998

e5

0.760

0.000

0.577

… ja jäägid

Kood 
residuals(mbf2, type="LDG2", verbose=FALSE, approx.z=!TRUE) |> (\(x) ("[<-"(x, upper.tri(x), NA)))() |> flexit()

X.

e1

e2

e3

e4

e5

e1

e2

0.007

e3

0.255

0.248

e4

0.658

1.791

0.616

e5

0.310

0.083

0.729

0.32

… ja mudeli sobitusastme indeksid

Kood 
data.frame(model=c("Algne", "Täiendatud"), rbind(fitses(mbf,3), fitses(mbf2, 3))) |> qflextable()

model

RMSEA

CFI

TLI

BIC

AIC

Algne

0.021

0.990

0.986

17 065.58

16 987.77

Täiendatud

0.011

0.997

0.996

17 044.41

16 960.12

St kui lubada 2. väite laadung spetsiifilisel faktoril, siis mudeli sobitusaste paraneb ning lokaalse sõltumatuse eeldust ei saa enam ümber lükata.

Kood 
anova(mbf, mbf2) |> flexit()

X.

AIC

SABIC

HQ

BIC

logLik

X2

df

p

mbf

16 987.77

17 027.45

17 015.09

17 065.58

-8 481.884

mbf2

16 960.12

17 003.10

16 989.71

17 044.41

-8 467.058

29.653

1

0

(See muudab natuke spetsiifilise faktori sisu: nüüd on liigsöömine vastandatud söömise piiramisele. Üldfaktori sisu muutub ka natuke: kõigi väidete kaalud muutuvad võrdsemaks.)

Invariantsus ja DIF

Cheung G. W., Rensvold R. B. (2002). Evaluating goodness-of-fit indexes for testing measurement invariance. Structural Equation Modeling: A Multidisciplinary Journal, 9(2), 233–255. Crossref. Web of Science.

Esineb sooline DIF. Siin sobib üsna hästi mudel, mis eeldab, et (1) 1.-3. väite puhul on meestel ja naistel sama nii eristusvõime kui ka raskusaste (2) 4. ja 5. väite raskusaste võivad olla erinevad

Kood 
data.frame(mudel=c("m2", "m1"), rbind( fitses(m1mgsj), fitses(m1mgs))) |> roundif(3) |> qflextable()

mudel

RMSEA

CFI

TLI

BIC

AIC

m2

0.024

0.949

0.944

17 003.74

16 919.45

m1

0.025

0.951

0.937

17 037.44

16 920.73

Kood 
 data.frame(mudel=c("m2", "m1"), anova(m1mgsj, m1mgs)) |> roundif(3) |> qflextable()

mudel

AIC

SABIC

HQ

BIC

logLik

X2

df

p

m2

16 919.45

16 962.43

16 949.05

17 003.74

-8 446.727

m1

16 920.73

16 980.24

16 961.71

17 037.44

-8 442.367

8.721

5

0.121

See tähendab, et “ankruteks” on 1.-3. väide ning 4. ja 5. väite puhul ühtlane DIF (erinev on ainult raskusaste).

Kood 
coefirt(m1mgsj) |> (\(x) data.frame(var=rownames(x), x))() |> roundif(3) |> qflextable()

var

X1.a

X1.b1

X1.b2

X2.a

X2.b1

X2.b2

e1

8.966

0.478

8.966

0.478

e2

1.523

1.423

1.523

1.423

e34

1.359

1.060

1.995

1.359

1.060

1.714

e5

1.892

3.126

1.892

2.964

3./4. väite iseloomulik kõver (NB! y-teljel on “oodatav skoor”, mitte iga vastusevariandi tõenäosus)

Kood 
itemplot(m1mgsj, item=3, type="score")

  1. väite iseloomulik kõver
Kood 
itemplot(m1mgsj, item=4)

Raskusastmete vahe: 4. väite puhul 0.271, 5. väite puhul 0.2.

(Tehniliselt mudelis on 3. ja 4. väide koos, kuid sellel väitel on kaks raskusastme parameetrit. St siin peab kohandama sõnastust, et mulje jääks õige - koefitsientide tabelis on 4 rida, aga väiteid on 5; joonise pealkirjas on “item 4”, aga see käib väidete tabeli viimase, viienda väite kohta)

Kood 
#https://pmc.ncbi.nlm.nih.gov/articles/PMC5505278/pdf/nihms855890.pdf
#https://search.r-project.org/CRAN/refmans/difR/html/difLogistic.html
#https://pmc.ncbi.nlm.nih.gov/articles/PMC2924271/

IRT mudelid

Ühefaktorilised mudelid: tabelis on eristusvõime parameeter a ja raskusastme parameeter b (klassikaline IRT parametriseering)

Kood 
cm1 <- coef(m1, IRTpar=TRUE, simplify=TRUE)
L2m1 <- mirt(eat2, itemtype="graded", verbose=FALSE)
cm2 <- coef(L2m1, IRTpar=TRUE, simplify=TRUE)
L3m1 <- mirt(eat3, itemtype="graded", verbose=FALSE)
cm3 <- coef(L3m1, IRTpar=TRUE, simplify=TRUE)
cm123 <- cbind(cm1$items[,1:3], " "=NA,  cm2$items[,1:3], "   " = NA, cm3$items[,1:3]) 
cm123 <- roundif(data.frame(" " = rownames(cm123), cm123),3)
names(cm123)<- c(" ", "L1_a", "L1_b", "L1_b2", "  ", "L2_a", "L2_b", "L2_b2", "   ", "L3_a", "L3_b", "L3_b2")

 qflextable(cm123 )

L1_a

L1_b

L1_b2

L2_a

L2_b

L2_b2

L3_a

L3_b

L3_b2

e1

10.189

0.454

14.405

0.489

11.366

0.506

e2

1.531

1.395

1.669

1.367

1.471

1.551

e34

1.374

1.029

1.765

1.388

1.163

1.836

1.310

1.096

1.851

e5

1.928

2.947

2.142

2.736

1.965

2.969

Bifaktori mudelite faktorlaadungid.

Kood 
ldngs <- function(x, pfx=NULL) {
  s <- summary(x, verbose=FALSE)
  res <- data.frame(s$rotF, s$h2)
  if(!is.null(pfx)) names(res) <- paste(pfx, names(res), sep="_")
  res
}
df <- data.frame("     "=names(eat1a), ldngs(mbf, "L1"), " "=NA, ldngs(L2mbf, "L2"), "  "=NA, ldngs(L3mbf, "L3"), check.names=FALSE)
In <- sapply(df, is.numeric)
df[, In] <- lapply(df[,In], function(x) ifelse(x==0, NA, x))
qflextable(roundif(df,3))

L1_G

L1_F1

L1_h2

L2_G

L2_F1

L2_h2

L3_G

L3_F1

L3_h2

e1

0.991

0.982

0.994

0.987

0.992

0.985

e2

0.671

0.451

0.704

0.496

0.656

0.430

e3

0.594

0.800

0.993

0.608

0.791

0.995

0.584

0.809

0.995

e4

0.685

0.726

0.996

0.675

0.735

0.996

0.667

0.742

0.996

e5

0.750

0.563

0.787

0.619

0.758

0.574

1F ja BF mudelite võrdlus

Korrelatsioonid ühefaktorilise mudeli ja BF mudeli üldfaktori skoori vahel (EAP = expected a posteriori).

Kood 
# * Vahemärkus. Miks on BF mudeli üldfaktori skoorid korreleeritud spetsiifiliste faktorite skooridega, on omaette küsimus. See on ilmselt expected a posteriori meetodi omapära. Mudelis on need faktorid üksteisest sõltumatud. Sõltumatud skoorid saame siis, kui kasutada meetodina "plausible" (sarnane nt PISA-testi plausible values teemaga), aga sel juhul tekivad mõned probleemid. Aga EAP meetod on skoorimise jaoks lihtsam (seda on ka sarnases valdkonnas kasutatud nt https://jpro.springeropen.com/articles/10.1186/s41687-022-00464-9), tõenäoliste väärtuste meetod annab iga kord erineva tulemuse (!) ja seetõttu tuleks iga analüüs läbi teha mitu korda (PISA soovituse järgi 10 korda, kuid vt ka https://pmc.ncbi.nlm.nih.gov/articles/PMC6425093/, mis ütleb, et miinimum on 20) ja raporteerida tulemuste keskmine. Tüütu. Teine huvitav asi on, et "plausible" skoorid on päriselt ka normaaljaotusega. # hist(fscores(mm2cbf, method="plausible")[,1]) 

# Plausible values kohta vt veel:
# 
# * https://link.springer.com/chapter/10.1007/978-3-030-47515-4_3
# 
# * https://ierinstitute.org/fileadmin/Documents/IERI_Monograph/Volume_2/IERI_Monograph_Volume_02_Chapter_01.pdf
# 
# * https://www.edmeasurementsurveys.com/IRT/Ch14.html

Võrdleme lihtsa ühefaktorilise mudeli (mudel1) ja esimese bifaktori mudeli (mudel 2 - spetsiifilisele faktorile laaduvad küsimused 3 ja 4) skoore.

Kood 
es1 <- fscores(m1, response.pattern=eat1, method="EAP")
es1a <- fscores(mbf, response.pattern=eat1a, method="EAP")
# ˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇˇ
es1b <- fscores(mbf2, response.pattern=eat1a, method="EAP")
# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
cor.test(es1[,1], es1a[,1])

    Pearson's product-moment correlation

data:  es1[, 1] and es1a[, 1]
t = 1619.5, df = 4907, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9990121 0.9991167
sample estimates:
      cor 
0.9990659
Kood 
plot(es1[,1], es1a[,1], xlab="Model 1 (1-factor)", ylab="Model 2 (bifactor)")
abline(lsfit(es1[,1], es1a[,1]), col="red")
text(-0.5, 2, "r=0.999")

Nüüd võrdleme esimese ja teise bifaktori mudeli skoore (mudel 2 vs. mudel 3 - esimeses on spetsiifilise faktoriga seotud küsimused 3 ja 4, teises 2,3 ja 4. )

Kood 
cor.test(es1a[,1], es1b[,1])

    Pearson's product-moment correlation

data:  es1a[, 1] and es1b[, 1]
t = 247.32, df = 4907, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9600147 0.9641730
sample estimates:
      cor 
0.9621498
Kood 
plot(es1a[,1], es1b[,1], xlab="Model 2 (bifactor)", ylab="Model 3 (bifactor+)")
abline(lsfit(es1a[,1], es1b[,1]), col="red")

Kood 
Li <- subset(L, incl)

Li$L1es <- es1[,1]
Li$L2es <- fscores(m1, response.pattern=eat2, method="EAP")[,1]
Li$L3es <- fscores(m1, response.pattern=eat3, method="EAP")[,1]
Kood 
op<-par(mfcol=c(1,3))
hist(Li$L1es, main="Wave 1")
hist(Li$L2es, main="Wave 2")
hist(Li$L3es, main="Wave 3")

Kood 
par(op)

Stabiilsus

Korrelatsioonid 3 laine vahel - diagonaali all ühefaktoriline mudel, diagonaalist üleval bifaktori mudeli üldfaktor.

Kood 
es.names <- c("L1es", "L2es", "L3es")
cor(Li[,es.names], use="pair") |> round(3) -> kors1
kors <- data.frame(" " = rownames(kors1), kors1, check.names = FALSE)
qflextable(kors)

L1es

L2es

L3es

L1es

1.000

0.654

0.602

L2es

0.654

1.000

0.642

L3es

0.602

0.642

1.000

Keskmiste erinevused…

Kood 
t.test(Li$L1es, Li$L3es, paired=TRUE)

    Paired t-test

data:  Li$L1es and Li$L3es
t = 2.6147, df = 2940, p-value = 0.008977
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.008348972 0.058418332
sample estimates:
mean difference 
     0.03338365
Kood 
t.test(Li$L2es, Li$L3es, paired=TRUE)

    Paired t-test

data:  Li$L2es and Li$L3es
t = 1.6722, df = 2376, p-value = 0.09461
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -0.003840642  0.048323705
sample estimates:
mean difference 
     0.02224153

St 3. laines on skoor oluliselt väiksem kui esimeses.

Kood 
t.test(Li$L1es, Li$L2es, paired=TRUE)

    Paired t-test

data:  Li$L1es and Li$L2es
t = 1.9785, df = 2548, p-value = 0.04798
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.0002265433 0.0507550722
sample estimates:
mean difference 
     0.02549081

KA esimese ja 2. laine vahel on väike erinevus.

Seosed KMI, vanuse ja sooga

KMI on arvutatud 1. ja 3. laine eneseraporteeritud kaalu ja pikkuse põhjal. Kaalu ja pikkuse eeltöötlus:

  • kaal alla 40 ja üle 200 kg asendatud NA-ga
  • pikkus alla 140 cm asendatud NA-ga
  • meetrites raporteeritud pikkus (3 juhul: 1.5 või 1.6) korrutatud 100-ga
  • KMI väärtused üle 50 ja alla 15 asendatud NA-ga

Korrelatsioonitabelis on (näidisena) ka KMI, soo ja vanuse keskmised korrelatsioonid 10 tõenäolise väärtusega (plausible value) nii ühefaktorilise (PV1F) kui ka bifaktori mudeli (PVBF) põhjal. Esimese hooga tundub mõstili eelistada EAP meetodit, mille puhul nn konvergentsed korrelatsioonid on suuremad.

Kood 
kaal <- Li$L1Kaal
# table(kaal)
kaal[kaal < 40] <- NA
kaal[kaal>200] <- NA
pikkus <- Li$L1Pikkus
#table(pikkus)
jupskid <- which(pikkus %in% c(1.5,1.6))
pikkus[jupskid] <- 100* pikkus[jupskid]
pikkus[pikkus<140] <- NA
pikkus <- pikkus/100

L1kmi <- kaal / pikkus^2
L1kmi[L1kmi > 50] <- NA
L1kmi[L1kmi < 15] <- NA

kaal <- Li$L3Kaal
# table(kaal)
kaal[kaal < 40] <- NA
kaal[kaal>200] <- NA
pikkus <- Li$L3Pikkus
# table(pikkus)
jupskid <- which(pikkus %in% c(1.5,1.6))
pikkus[jupskid] <- 100* pikkus[jupskid]
pikkus[pikkus<140] <- NA
pikkus <- pikkus/100

L3kmi <- kaal / pikkus^2
L3kmi[L3kmi > 50] <- NA
L3kmi[L3kmi < 15] <- NA

Li$L1kmi <- L1kmi
Li$L3kmi <- L3kmi

taust <- data.frame(KMI = L1kmi, Sugu = Li$RR.Sugu, Vanus = as.numeric(Li$RR.Vanus))
k1 <- cor(taust, Li[, c("L1es")], use="com")
taust3 <- data.frame(KMI = L3kmi, Sugu = Li$RR.Sugu, Vanus = as.numeric(Li$RR.Vanus))
k3 <- cor(taust3, Li[, c("L3es")], use="com")

df <- data.frame(" "=rownames(k1), k1, "  "=NA, k3)  |> roundif(3)
names(df) <- c(" ", "L1", "      ", "L3")
qflextable(df)

L1

L3

KMI

0.236

0.241

Sugu

0.110

0.096

Vanus

-0.218

-0.233

Seosed isiksuse, EEK-2 ja impulsiivsusega

  • Isiksuse skaala: 20 väitega lühiversioon 100NP põhjal

  • Impulsiivsus: negative & positive urgency, Short version of the UPPS-P Impulsive Behavior Scale (SUPPS-P; Cyders, Littlefield, Coffey, & Karyadi, 2014); SUPPS-P (Cyders et al., 2014)

  • NB! Isiksuse andmed 3. lainest, impulsiivsus 2. lainest. Kõik korrelatsioonid on olulised (p< 0.001)

Kood 
NEOAC2 <- Li[, paste0("L3ISIKSUS", 1:22)] #  grab(Li, SRC %in% "100NP")#[c("N", "E", "O", "A", "C")]
key(NEOAC2) <- key(L, SRC %in% "100NP")
NEOAC2$L3ISIKSUS6 <- 7-NEOAC2$L3ISIKSUS6
NEOAC2 <- compute_scores(NEOAC2, SCALE %in% c("N", "E", "O", "A", "C"))[, c("N", "E", "O", "A", "C")]
NEOACS <- apply(NEOAC2, 2, scale)
IMP <- Li[, paste0("L2C", 35:42)] # SRC %in% "IMP")
IMP[] <- lapply(IMP, \(x) 5-x)
pu <- rowSums(IMP[,1:4])
nu <- rowSums(IMP[,5:8])

imps <- sub("[", "", sapply(strsplit(key(L, SRC %in% "IMP")$LAB, "]"), "[[", 1), fixed=TRUE)

qflextable(data.frame(nr=1:8, Väide=imps))

nr

Väide

1

Väga heas tujus olles kipun sattuma olukordadesse, mis võivad mulle probleeme tekitada.

2

Väga heas tujus olles kipun kaotama kontrolli.

3

Kui olen väga elevil, teen asju, mis panevad teisi muretsema või ehmatavad neid.

4

Kui olen tõeliselt elevil, kipun tegutsema mõtlematult.

5

Kui tunnen end halvasti, teen hetke ajel enesetunde parandamiseks sageli asju, mida hiljem kahetsen.

6

Vahel, kui tunnen end halvasti, ei suuda ma oma tegevust lõpetada, kuigi see teeb mu enesetunde halvemaks.

7

Kui mu meeleolu on halb, tegutsen sageli mõtlematult.

8

Kui tunnen end tõrjutuna, ütlen sageli asju, mida hiljem kahetsen.

Kood 
IMPS <- data.frame( PU=pu, NU=nu )
Kood 
cor_neo <- psych::corr.test(NEOACS, Li[, c("L3es")], use="pair")
cor_imp <- psych::corr.test(IMPS, Li[, c("L2es")], use="pair")
bindit <- function(x) setNames(data.frame(x$r , N=x$n[,1]), c("eat_1F", "N"))
# 
rbind(bindit(cor_neo), bindit(cor_imp)) |> flexit()

X.

eat_1F

N

N

0.283

3 872

E

-0.049

3 870

O

-0.064

3 886

A

-0.112

3 891

C

-0.063

3 878

PU

0.165

3 100

NU

0.241

3 074

Ülesöömise skoori (1F mudel) prediktorid lineaarse mudeli põhjal:

Kood 
lms <- lm(Li$L1es ~ as.matrix(NEOAC2) + Li$RR.Sugu + Li$RR.Vanus) |> summary()
df <- data.frame(x = c("Vabaliige", "N", "E", "O", "A", "C", "Sugu (1=M,2=N)", "Vanus"), coef(lms))
names(df) <- c("  ", "b", "SE(b)", "t", "p")
df[,2:4] <- round(df[,2:4],3)
df$p<-round(df$p,5)
qflextable(df)

b

SE(b)

t

p

Vabaliige

-0.044

0.143

-0.309

0.75762

N

0.044

0.004

10.170

0.00000

E

0.006

0.006

1.127

0.25965

O

-0.001

0.005

-0.201

0.84049

A

-0.005

0.007

-0.648

0.51699

C

-0.013

0.005

-2.403

0.01635

Sugu (1=M,2=N)

0.120

0.031

3.811

0.00014

Vanus

-0.008

0.001

-9.310

0.00000

(R2 = 0.1057009)

EEK-i skooride (depressioon, üldistunud ärevushäire, asteenia, unehäired) korrelatsioonid ülesöömise skooridega 3 laines:

Kood 
eeks <- c("EEK.dep", "EEK.anx", "EEK.ast", "EEK.ins")
cor_eek <- data.frame(L1=cor(Li[, paste0("L1", eeks)], Li[, c("L1es")], use="com"), 
L2=cor(Li[, paste0("L2", eeks)], Li[, c("L2es")], use="com"), 
L3=cor(Li[, paste0("L3", eeks)], Li[, c("L3es")], use="com"), check.names=FALSE)
rownames(cor_eek) <- eeks
cor_eek |> flexit()

X.

L1

L2

L3

EEK.dep

0.278

0.305

0.249

EEK.anx

0.307

0.328

0.286

EEK.ast

0.266

0.291

0.246

EEK.ins

0.174

0.207

0.171

Diagnooside ennustamine

Järgnevad analüüsid on raporteeritud funktsiooniga jtools::summ: log.regressiooni kordajad on tabelis nn “exp(x)” kujul, st OR. Prediktorid (v.a sugu) on standardiseeritud, sh vanus ja vanuse ^2.

Kood 
Ri<-subset(R, incl)
Li$F32 <- (Ri$F32_6 > 0)*1
Li$F33 <- (Ri$F33F34_6 > 0)*1
Li$dep <- Li$F32 | Li$F33
Li$F41 <- (Ri$F41_6 > 0)
Li$F32a <- (Ri$F32_2021>0)*1
Li$F33a <- (Ri$F33F34_2021>0)*1
Li$depa <- Li$F32a | Li$F33a
Li$F41a <- (Ri$F41_2021>0)*1

g10 <- glm(dep ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1es, data=Li, family=binomial) 
g11 <- glm(dep ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1EEK.dep + L1es, data=Li, family=binomial) 
g12 <- glm(dep ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1EEK.dep + L1EEK.ast + L1es, data=Li, family=binomial) 
g20 <- glm(depa ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1es , data=Li, family=binomial)
g21 <- glm(depa ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1EEK.dep + L1es , data=Li, family=binomial) 
g22 <- glm(depa ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1EEK.dep + L1EEK.ast + L1es , data=Li, family=binomial) 
g30 <- glm(F41 ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1es, data=Li, family=binomial) 
g31 <- glm(F41 ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1EEK.anx + L1es, data=Li, family=binomial) 
g40 <- glm(F41a ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1es, data=Li, family=binomial)
g41 <- glm(F41a ~ RR.Sugu + RR.Vanus + I(RR.Vanus^2) + L1EEK.anx + L1es, data=Li, family=binomial)

Sõltuv muutuja: depressiooni (F32-34) 5 aasta levimus (diagn 2016-20)

1) Söömishäire riskiskoor on oluline prediktor

Kood 
summ1(g10)
χ²(4) 100.6170
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.0400
Pseudo-R² (McFadden) 0.0289
AIC 3386.1106
BIC 3418.6048
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0821 0.0693 0.0972 -28.9935 0.0000
RR.Sugu 1.8399 1.5059 2.2480 5.9652 0.0000
RR.Vanus 1.6263 0.9894 2.6734 1.9178 0.0551
I(RR.Vanus^2) 0.6614 0.4029 1.0858 -1.6344 0.1022
L1es 1.3620 1.2490 1.4851 6.9961 0.0000
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

2) Söömishäire riskiskooril on lisavaliidsus (incremental validity)

Kood 
summ1(g11)
χ²(5) 393.4799
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.1596
Pseudo-R² (McFadden) 0.1193
AIC 2916.7579
BIC 2955.4424
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0659 0.0546 0.0795 -28.4368 0.0000
RR.Sugu 1.8375 1.4831 2.2766 5.5648 0.0000
RR.Vanus 3.7330 2.1946 6.3498 4.8601 0.0000
I(RR.Vanus^2) 0.3317 0.1958 0.5620 -4.1022 0.0000
L1EEK.dep 2.1780 1.9886 2.3855 16.7678 0.0000
L1es 1.1514 1.0461 1.2672 2.8811 0.0040
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

3) … lisavaliidsus v.2.0 (kontrollime lisaks EEK-2 depressiooniskoorile ka asteenia skooir)

Kood 
summ1(g12)
χ²(6) 411.0578
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.1692
Pseudo-R² (McFadden) 0.1268
AIC 2844.6647
BIC 2889.6740
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0656 0.0543 0.0793 -28.1312 0.0000
RR.Sugu 1.7739 1.4283 2.2031 5.1846 0.0000
RR.Vanus 4.1922 2.4482 7.1783 5.2228 0.0000
I(RR.Vanus^2) 0.2966 0.1740 0.5056 -4.4663 0.0000
L1EEK.dep 1.8335 1.6207 2.0743 9.6290 0.0000
L1EEK.ast 1.3310 1.1662 1.5192 4.2388 0.0000
L1es 1.1314 1.0261 1.2475 2.4759 0.0133
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

Sõltuv muutuja: depressiooni (F32-34) 2 kuu levimus (diagnoos 2021 (aanuar-veebruar)

1) Söömishäire riskiskoor on oluline prediktor

Kood 
summ1(g20)
χ²(4) 24.2907
p 0.0001
Pseudo-R² (Cragg-Uhler) 0.0258
Pseudo-R² (McFadden) 0.0232
AIC 1034.4219
BIC 1066.6493
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0142 0.0096 0.0211 -21.2665 0.0000
RR.Sugu 2.0108 1.2852 3.1460 3.0585 0.0022
RR.Vanus 2.4311 0.8107 7.2904 1.5853 0.1129
I(RR.Vanus^2) 0.3870 0.1266 1.1833 -1.6648 0.0959
L1es 1.3064 1.0893 1.5668 2.8819 0.0040
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

2) Söömishäire riskiskooril ei ole lisavaliidsust

Kood 
summ1(g22)
χ²(6) 132.3044
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.1494
Pseudo-R² (McFadden) 0.1358
AIC 855.7129
BIC 900.3375
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0092 0.0058 0.0144 -20.2927 0.0000
RR.Sugu 1.9707 1.2219 3.1783 2.7819 0.0054
RR.Vanus 7.8569 2.4508 25.1878 3.4681 0.0005
I(RR.Vanus^2) 0.1378 0.0420 0.4523 -3.2687 0.0011
L1EEK.dep 2.3451 1.8390 2.9904 6.8720 0.0000
L1EEK.ast 1.1085 0.8400 1.4629 0.7278 0.4667
L1es 1.0335 0.8468 1.2613 0.3241 0.7459
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

Sõltuv muutuja: “muude ärevushäirete” 5 aasta levimus (2016-2020)

1) Söömishäire riskiskoor on oluline prediktor

Kood 
summ1(g30)
χ²(4) 55.6488
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.0257
Pseudo-R² (McFadden) 0.0196
AIC 2789.2496
BIC 2821.7437
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0597 0.0491 0.0725 -28.4094 0.0000
RR.Sugu 1.8429 1.4637 2.3203 5.2015 0.0000
RR.Vanus 0.8264 0.4821 1.4165 -0.6936 0.4879
I(RR.Vanus^2) 1.3138 0.7691 2.2442 0.9990 0.3178
L1es 1.2393 1.1217 1.3692 4.2182 0.0000
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

2) Söömishäire riskiskooril ei ole lisavaliidsust

Kood 
summ1(g31)
χ²(5) 208.6424
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.0990
Pseudo-R² (McFadden) 0.0770
AIC 2514.2852
BIC 2553.0426
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0542 0.0441 0.0666 -27.6198 0.0000
RR.Sugu 1.6173 1.2702 2.0592 3.9004 0.0001
RR.Vanus 1.1481 0.6590 2.0001 0.4875 0.6259
I(RR.Vanus^2) 1.0617 0.6119 1.8421 0.2128 0.8315
L1EEK.anx 1.9331 1.7414 2.1459 12.3707 0.0000
L1es 1.0419 0.9347 1.1615 0.7409 0.4587
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

Sõltuv muutuja: “Muude ärevushäirete” 2 kuu levimus (2021, jaanuar-veebruar)

1) Söömishäire riskiskoor ei ole oluline prediktor

Kood 
summ1(g40)
χ²(4) 15.6013
p 0.0036
Pseudo-R² (Cragg-Uhler) 0.0201
Pseudo-R² (McFadden) 0.0184
AIC 843.3180
BIC 875.5454
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0121 0.0079 0.0185 -20.3424 0.0000
RR.Sugu 1.7442 1.0631 2.8615 2.2023 0.0276
RR.Vanus 0.5347 0.1789 1.5987 -1.1203 0.2626
I(RR.Vanus^2) 1.5669 0.5015 4.8963 0.7726 0.4398
L1es 1.2214 0.9910 1.5053 1.8753 0.0608
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.

2) Söömishäire riskiskooril ei ole lisavaliidsust

Kood 
summ1(g41)
χ²(5) 75.8337
p 0.0000
Pseudo-R² (Cragg-Uhler) 0.0997
Pseudo-R² (McFadden) 0.0918
AIC 762.3418
BIC 800.7727
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.0103 0.0065 0.0161 -19.8304 0.0000
RR.Sugu 1.3892 0.8380 2.3032 1.2746 0.2025
RR.Vanus 0.9272 0.3056 2.8132 -0.1335 0.8938
I(RR.Vanus^2) 1.0407 0.3265 3.3174 0.0675 0.9462
L1EEK.anx 2.2898 1.8623 2.8156 7.8569 0.0000
L1es 0.9887 0.7917 1.2348 -0.0999 0.9204
Standard errors: MLE; Continuous predictors are mean-centered and scaled by 1 s.d. The outcome variable remains in its original units.
Kood 
save(Li, file="Li.rda")