Re: [R] cross-classified random factors in lme without blocking - success
Sundar's advice seems to do the trick. Here is a small simulation of a cross-classified random model with extraction of the fitted variance components from lme: library(nlme) set.seed(18112003) na - 20 nb - 20 sigma.a - 2 sigma.b - 3 sigma.res - 1 mu - 0.5 n - na*nb a - gl(na,1,n) b - gl(nb,na,n) u - gl(1,1,n) ya - rnorm(na,mean=0,sd=sigma.a) yb - rnorm(nb,mean=0,sd=sigma.b) y - model.matrix(~a-1) %*% ya + model.matrix(~b-1) %*% yb + rnorm(n,mean=mu,sd=sigma.res) fm - lme(y~1,random=list(u=pdBlocked(list(pdIdent(~a-1),pdIdent(~b-1) v - sqrt(diag(as.matrix(fm$modelStruct$reStruct$u))[c(1,na+1)]) v - fm$sigma * c(v, 1) names(v) - c(Sigma.a,Sigma.b,Sigma.res) print(v) Sigma.a Sigma.b Sigma.res 2.450700 3.650427 1.046689 Gordon At 05:29 AM 16/11/2003, Sundar Dorai-Raj wrote: Gordon Smyth wrote: Many thanks for help from Peter Dalgaard, Douglas Bates and Bill Venables. As a result of their help, here is a working example of using lme to fit an additive random effects model. The model here is effectively y~a+b with a and b random: y - rnorm(12) a - gl(4,1,12) b - gl(3,4,12) u - gl(1,1,12) library(nlme) fm - lme(y~1,random=list(u=pdBlocked(list(pdIdent(~a-1),pdIdent(~b-1) summary(fm) I still can't see though how to extract the three variance components (for a and b and for the residual) from the fitted object 'fm'. Thanks Gordon Would this work for you? It still needs some parsing to extract sigma_a and sigma_b, but that should be simple. v - as.matrix(fm$modelStruct$reStruct$u) c(sqrt(diag(v)), Residual = 1) * fm$sigma Regards, Sundar __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] cross-classified random factors in lme without blocking
Many thanks for help from Peter Dalgaard, Douglas Bates and Bill Venables. As a result of their help, here is a working example of using lme to fit an additive random effects model. The model here is effectively y~a+b with a and b random: y - rnorm(12) a - gl(4,1,12) b - gl(3,4,12) u - gl(1,1,12) library(nlme) fm - lme(y~1,random=list(u=pdBlocked(list(pdIdent(~a-1),pdIdent(~b-1) summary(fm) I still can't see though how to extract the three variance components (for a and b and for the residual) from the fitted object 'fm'. Thanks Gordon __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] cross-classified random factors in lme without blocking
Gordon Smyth [EMAIL PROTECTED] writes: On page 165 of Mixed-Effects Models in S and S-Plus by Pinheiro and Bates there is an example of using lme() in the nlme package to fit a model with crossed random factors. The example assumes though that the data is grouped. Is it possible to use lme() to fit crossed random factors when the data is not grouped? E.g., y - rnorm(12); a=gl(4,1,12); b=gl(3,4,12). Can I fit an additive model like y~a+b but with a and b random? Everything I've tried gives an error: lme(y~1,random=~1|(a+b)) Error in switch(mode(object), name = , numeric = , call = object, character = as.name(object), : [[ cannot be of mode ( lme(y~1,random=~a+b) Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups A standard trick is to define a grouping with one level: one - rep(1,length(y) lme(, random=~pdBlocked(.)|one) (Sorry, I'm a little rusty on the syntax, but just follow the example in PB) AFAIR, it also works with random=list(a=~1,one=~b) and vice versa. (The model is the same but you get different DF calculations, none of which are correct in the completely balanced case...) -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] cross-classified random factors in lme without blocking
Peter Dalgaard [EMAIL PROTECTED] writes: Gordon Smyth [EMAIL PROTECTED] writes: On page 165 of Mixed-Effects Models in S and S-Plus by Pinheiro and Bates there is an example of using lme() in the nlme package to fit a model with crossed random factors. The example assumes though that the data is grouped. Is it possible to use lme() to fit crossed random factors when the data is not grouped? E.g., y - rnorm(12); a=gl(4,1,12); b=gl(3,4,12). Can I fit an additive model like y~a+b but with a and b random? Everything I've tried gives an error: lme(y~1,random=~1|(a+b)) Error in switch(mode(object), name = , numeric = , call = object, character = as.name(object), : [[ cannot be of mode ( lme(y~1,random=~a+b) Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups A standard trick is to define a grouping with one level: one - rep(1,length(y) lme(, random=~pdBlocked(.)|one) (Sorry, I'm a little rusty on the syntax, but just follow the example in PB) AFAIR, it also works with random=list(a=~1,one=~b) and vice versa. Not sure about that. (The model is the same but you get different DF calculations, none of which are correct in the completely balanced case...) I realize that it is awkward to use lme to fit models with crossed random effects. As Saikat DebRoy and I described in a recent preprint http://www.stat.wisc.edu/~bates/reports/MultiComp.pdf we now have a good handle on the computational methods for mixed-effects models with nested or crossed or partially crossed random effects. Both the nlme and the lme4 packages are based on structures that are tuned to nested random effects and do not easily accomodate crossed random effects. I have a draft of the contents of classes and methods for fitting linear mixed-effects models with nested or crossed or ... but it is a long way from the draft to working, tested code. Although it will take some time to get all the pieces in place I do offer some encouragement that this awkward phrasing of crossed random effects will some day be behind us. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] cross-classified random factors in lme without blocking
Douglas Bates [EMAIL PROTECTED] writes: (Sorry, I'm a little rusty on the syntax, but just follow the example in PB) AFAIR, it also works with random=list(a=~1,one=~b) and vice versa. Not sure about that. Sorry. It's certainly not correct as written. It has to be something like list(a=1,one=pdIdent(form=~b-1)) otherwise you get a general symmetric covariance for the effect of b. (The model is the same but you get different DF calculations, none of which are correct in the completely balanced case...) I realize that it is awkward to use lme to fit models with crossed random effects. As Saikat DebRoy and I described in a recent preprint http://www.stat.wisc.edu/~bates/reports/MultiComp.pdf .../MixedComp.pdf, right? we now have a good handle on the computational methods for mixed-effects models with nested or crossed or partially crossed random effects. Both the nlme and the lme4 packages are based on structures that are tuned to nested random effects and do not easily accomodate crossed random effects. I have a draft of the contents of classes and methods for fitting linear mixed-effects models with nested or crossed or ... but it is a long way from the draft to working, tested code. Although it will take some time to get all the pieces in place I do offer some encouragement that this awkward phrasing of crossed random effects will some day be behind us. Looking forward to it... :-) -- O__ Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
Re: [R] cross-classified random factors in lme without blocking
Douglas Bates [EMAIL PROTECTED] writes: ... I realize that it is awkward to use lme to fit models with crossed random effects. As Saikat DebRoy and I described in a recent preprint http://www.stat.wisc.edu/~bates/reports/MultiComp.pdf we now have a good handle on the computational methods for mixed-effects models with nested or crossed or partially crossed random effects. That URL should have been http://www.stat.wisc.edu/~bates/reports/MixedComp.pdf __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
