God ide. Måske skulle vi lave et forum af en art i stedet. Måske findes der noget allerede (på dansk). Hvad har du i tankerne om det konkrete indhold?

/Kristian

Jeg foreslår, at vi bare bruger forummet her. Det gør det selvfølgelig svært at sende screendumps rundt, men måske kunne vi reservere en mailingliste til det.

/Kristian

Hej alle,

Så løb jeg ind i det første seriøse problem i R. Jeg har estimeret en lineær multilevel model. Vi ser her på individ [i]i[/i] til tiden [i]t[/i]. Jeg ser på indkomstudviklingen over livet.

INDKOMST = INTERCEPT + ALDER + ÅRSLEDIGHEDSGRAD

På multilevelform er modellen:

Niveau 1
yit = b0i + b1ix1 + b2ix2 + eit

Niveau 2
b0i = gamma00 + u0i
b1i = gamma10 + u1i
b2i = gamma20 + u2i

Her antages det at uki, hvor k = {0,1,2}, er multivariat normalfordelt, dvs. trivariat fordelt.

Koden for modellen i R, er:

[i][b]m1 m1
Linear mixed-effects model fit by maximum likelihood
Data: dstclean
Log-likelihood: -14766.62
Fixed: INDK ~ ALDER + ARL
(Intercept) ALDER ARL
11.79998411 0.02533442 -0.18167027

Random effects:
Formula: ~1 + ALDER + ARL | LBNR
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.82344428 (Intr) ALDER
ALDER 0.04542280 -0.696
ARL 0.50494872 -0.371 0.073
Residual 0.28588168

Number of Observations: 34566
Number of Groups: 3194 [/i]

Mit problem er, at jeg ikke kan finde ud af at trække varianserne og kovarianserne ud af matricen for random komponenterne på niveau 2 (jf. den trivariate n-fordeling ovenfor).

Jeg har nemlig brug for at kunne trække alle elementer ud af matricen, da jeg skal bruge dem til at regne videre senere hen.

Jeg har prøvet med:
[i]
[b]m1$modelStruct$corStruct[1][/b][/i]

Har også prøvet at lave summary() på m1, men her kan jeg heller ikke finde ud af trække noget ud. str(m1) for modelStruct ser sådan ud:

[i]> str(m1$modelStruct)
List of 1
$ reStruct:List of 1
..$ LBNR:Classes `pdLogChol`, `pdSymm`, `pdMat` atomic [1:6] 1.058 -2.171 0.455 -0.111 -0.655 ...
.. .. ..- attr(*, "formula")=Class `formula` length 2 ~1 + ALDER + ARL
.. .. .. .. ..- attr(*, ".Environment")=length 2
.. .. ..- attr(*, "Dimnames")=List of 2
.. .. .. ..$ : chr [1:3] "(Intercept)" "ALDER" "ARL"
.. .. .. ..$ : chr [1:3] "(Intercept)" "ALDER" "ARL"
..- attr(*, "settings")= num [1:4] 0 1 0 4
..- attr(*, "class")= chr "reStruct"
..- attr(*, "plen")= Named int 6
.. ..- attr(*, "names")= chr "LBNR"
- attr(*, "settings")= num [1:4] 0 0 1 4
- attr(*, "class")= chr [1:3] "lmeStructInt" "lmeStruct" "modelStruct"
- attr(*, "pmap")= logi [1:6, 1] TRUE TRUE TRUE TRUE TRUE TRUE
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr "reStruct"
> [/i]

Hvis der er nogen, der kan lure, hvordan jeg får mine estimater for varianser og kovarianser trukket ud, så er jeg mere end lykkelig. Har prøvet at kigge på r-help, men kan ikke rigtig finde noget. Jeg er ikke så habil ud i R, så det kan også være, at jeg ganske enkelt overser noget.

vh Kristian

Koden for min lme er:

m1

undskyld? den fangede jeg ikke helt ;)

Jamen så nåede jeg ikke at svare.

Laver du en grafisk fremstilling samtidigt i R?

Hvilket datasæt anvender du - jeg kan se at det er et enorm stort antal respondendter.

Mange hilsener
Christian Christrup Kjeldsen