[Sorry, In case this is repeating a message already sent to the list].
I am trying to move a project to R (base or nlme), for which I have a
partial
solution in BMDP 8V. Here is the 8V control language:
/input title='Augenbewegungen'.
variables=4.
file='latm.dat'.
format='11x,f7.0,f7.0,f12.4,f10.0'.
/variables namesllat,rlat,vg,diff.
/design
ndep = 4.
level=21,2,2,9,3.
names = person,retest,direct,flicker,ground.
random = person.
fixed = retest,ground,direct,flicker.
model = 'P,R,G,D,F'.
/end.
/finish.
BMDP 8V calculates the anova and reports F-values for tests of all fixed
effect
interactions against the error from the interaction with person
variables.
Also estimates of variance components. (BMDP is used on uVAX3100 , 24MB,
VMS system, CPU time 20 s for this problem). I tried to translate the
setup
to R. The nearest equivalent I got was:
latm <- read.table("latm.dat" )
latm.aov <- aov( llat ~ person + retest * direct * flicker * ground
+ Error ( person : ( retest * direct * flicker * ground )))
summary(latm.aov) at least reports the same F-values as 8V does. Though
the
R computations thrash my system (FreeBSD, K6-233, 64MB). I have to
triple the swapspace before starting, and I wait about an hour for
either rebooting the system or having a lucky outcome. There seems to be
something fundamentally wrong.
It seems I have to move to nmle to get variance components. I did not
manage to translate my problem to mle language yet.
Thanks for any help
--
Dipl.-Math. Wilhelm Bernhard Kloke
Institut fuer Arbeitsphysiologie an der Universitaet Dortmund
Ardeystrasse 67, D-44139 Dortmund, Tel. 0231-1084-257
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
"Wilhelm B. Kloke" <wb at vestein.arb-phys.uni-dortmund.de> writes:> latm <- read.table("latm.dat" ) > latm.aov <- aov( llat ~ person + retest * direct * flicker * ground > + Error ( person : ( retest * direct * flicker * ground ))) > > summary(latm.aov) at least reports the same F-values as 8V does. Though > the > R computations thrash my system (FreeBSD, K6-233, 64MB). I have to > triple the swapspace before starting, and I wait about an hour for > either rebooting the system or having a lucky outcome. There seems to be > something fundamentally wrong. > > It seems I have to move to nmle to get variance components. I did not > manage to translate my problem to mle language yet.(That's nlme/lme - (non-)linear mixed effect) aov() does take rather a long time on problems of this size (2268 observations, 16 error strata). I tried it with simulated data and saw the the memory footprint go to about 300MB for a while... The algorithm in R is not as good at taking advantage of balancedness between error terms as programs like Genstat and does the orthogonal decomposition by brute force, leading to linear problems with a huge number of predictors. Of course, given sufficient knowledge of the particular model, you can optimize very heavily and everything reduces to means and sums of squares thereof. That is probably what BMDP 8V does. Variance components can be obtained from the estimated stratum error terms by inverting the linear relation between them, which is a bit tedious to work out, but I suspect you'll find it even more tedious to set up the model in lme()... The latter involves explicitly setting up the block diagonal covariance matrix of random effects, with each block proportional to the identity matrix. Something like random=list(person=pdBlocked( list(pdIdent(~retest-1), ...., pdIdent(~retest:direct:flicker:ground-1)))) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._