Dear John, R-sig-mixed-models is more suited for this kind of questions. All follow-up mail should be posted only to that mailing list.
It seems like varIdent() by default relevels the grouping factor and that the user cannot control this. Best regards, ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium + 32 2 525 02 51 + 32 54 43 61 85 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -----Oorspronkelijk bericht----- Van: R-help [mailto:r-help-boun...@r-project.org] Namens Kornak, John Verzonden: dinsdag 23 december 2014 22:29 Aan: r-help@R-project.org Onderwerp: [R] nlme package: changing reference for varIdent parameter estimates in summary.gls Dear R experts, I am running gls models with heterogeneous group variances using the glm function in the nlme package with varIdent weights. I am interested in controlling the baseline level for the group variance function standard deviation estimates by applying the relevel function to the group variable. When I use relevel, the baseline level in the coefficient table changes as expected, but the baseline level does not change for the variance function table; for example, see the fitted models gls1 vs. gls2 in the contrived example below. Does anyone have a suggestion as to how I can change the baseline level for the variance function output please? In addition to the example below, I have tried specifying the value argument as per the varIdent help page, e.g. variIdent(c(no=0.5), form = ~1|group) and have google searched / checked help pages for solutions without success. I am running R version 3.1.0 on an iMac OSX v. 10.9.5 Thank you in advance John Kornak > library(nlme) > group <- factor(c(rep("no",20),rep("yes",20))) > set.seed(2) > outcome <- c(rnorm(20,0,2),rnorm(20,5,4)) dataTest <- > data.frame(outcome,group) # Original model fit before releveling > gls1 <- gls(outcome ~ group, weights=varIdent(form = ~1|group), > data=dataTest) > summary(gls1) snip Variance function: Structure: Different standard deviations per stratum Formula: ~1 | group Parameter estimates: no yes 1.00000 2.23034 Coefficients: Value Std.Error t-value p-value (Intercept) 0.390922 0.4734001 0.825775 0.4141 groupyes 4.607951 1.1571140 3.982279 0.0003 snip Residual standard error: 2.11711 Degrees of freedom: 40 total; 38 residual # relevel the group so that yes is the reference > dataTest$group <- relevel(dataTest$group,"yes") > gls2 <- gls(outcome ~ group, weights=varIdent(form = ~1|group), > data=dataTest) > summary(gls2) snip Variance function: Structure: Different standard deviations per stratum Formula: ~1 | group Parameter estimates: no yes 1.00000 2.23034 ### no" is still the reference group here for the variance function Coefficients: Value Std.Error t-value p-value (Intercept) 4.998873 1.055843 4.734484 0e+00 groupno -4.607951 1.157114 -3.982279 3e-04 # yes has become the reference for the coefficients snip Residual standard error: 2.11711 Degrees of freedom: 40 total; 38 residual > --------------------------------------------------- John Kornak, PhD Associate Professor in Residence Department of Epidemiology and Biostatistics University of California, San Francisco Mission Hall: Global Health & Clinical Sciences Building 550 16th St, 2nd floor, Box #0560 San Francisco, CA 94158-2549 Tel: 415-514-8028 Fax: 415-514-8150 Email: john.kor...@ucsf.edu<mailto:john.kor...@ucsf.edu> [[alternative HTML version deleted]] Disclaimer Bezoek onze website / Visit our website<https://drupal.inbo.be/nl/disclaimer-mailberichten-van-het-inbo> ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
