Application of coxme
Hi David, coxme now uses the lme4 syntax or you have to specify fixed and random: coxme (Surv(time,status) ~ group + (1|nest), data=females) # or coxme (fixed=Surv(time,status) ~ group, random=~1|nest, data=females) Jonas
On Wed 05 Jun 2013 12:16:18 PM CEST, davidcostantini at libero.it wrote:
Dear All I am trying to use the coxme function to compare survival among 6 experimental groups. I have a fixed factor (group), a random factor (brood), and a censoring variable (0 = still alive; 1 = dead). My response variable is the age of the dead individual or of the individual still alive at the end of the experiment. I have tried to use this function: coxme (Surv(time,status) ~ group, random=~1|nest, data=females) Time would be the age, while status is the censoring variable. R gives me an output and says that the random argument of coxme is depreciated. I wonder if anyone of you may let me know if this function is correct and why the random factor is depreciated. Cheers David
----Messaggio originale---- Da: r-sig-mixed-models-request at r-project.org Data: 05/06/2013 12.00 A: <r-sig-mixed-models at r-project.org> Ogg: R-sig-mixed-models Digest, Vol 78, Issue 11 Send R-sig-mixed-models mailing list submissions to r-sig-mixed-models at r-project.org To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models or, via email, send a message with subject or body 'help' to r-sig-mixed-models-request at r-project.org You can reach the person managing the list at r-sig-mixed-models-owner at r-project.org When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-mixed-models digest..." Today's Topics: 1. Re: nlme: varIdent (Ben Bolker) 2. Re: Autocorrelation in a GAMM for nightly time series (Ben Bolker) 3. Re: Negative Variance (John Maindonald) ---------------------------------------------------------------------- Message: 1 Date: Wed, 5 Jun 2013 01:14:37 +0000 (UTC) From: Ben Bolker <bbolker at gmail.com> To: r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] nlme: varIdent Message-ID: <loom.20130605T031041-601 at post.gmane.org> Content-Type: text/plain; charset=us-ascii D Gill <oceancurrents at ...> writes: Hello, I am a student working on fisheries data using a mixed linear model. I have a question about the varIdent function. I have information on 140 fishers at nine sites, where the only random effect is for the site. As I do not have equal variances between sites, is it correct to use the varIdent argument for the sites? Do I lose any information by doing it this way? M1<-lme(resp~......,random = ~1|Site,weights = varIdent(form =~1|Site) )) This seems reasonable to me. However, note that you might want to consider (variable*site) interactions of any predictor variables that vary within sites (see e.g. Schielzeth and Forstmeier 2009) [e.g. if your measurement is total catch, effort is a covariate, and CPUE varies among sites, you would want effort|Site in your random effects specification]. I'm assuming you have a single measurement for each fisher -- otherwise you should also add a random effect for fishers ... Ben Bolker ------------------------------ Message: 2 Date: Wed, 5 Jun 2013 01:35:51 +0000 (UTC) From: Ben Bolker <bbolker at gmail.com> To: r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] Autocorrelation in a GAMM for nightly time series Message-ID: <loom.20130605T032637-397 at post.gmane.org> Content-Type: text/plain; charset=us-ascii Andrew Digby <andrewdigby at ...> writes: I'm analysing call counts of a nocturnal species in response to temporal and environmental effects. Over a period of 3 years I divide each night into 10 equal-length blocks, and count the calls of each sex in each block. [snip] gam(ncalls ~ sex + year + s(month) + s(time of night) + s(moon) + s(temperature) + s(rain) + offset(blocklength), family=negbin(c(1,5)), data=dc) Where 'time of night' = 1-10 = the block number during the night, and blocklength = the length of that block (this changes throughout the year; hence the offset). Not surprisingly, the acf shows significant correlation between adjacent time blocks, with a periodicity of 10. This is because counts in a particular block will be similar to those in the blocks immediately before and after, and also to the same block on other nights (the species has a quite regular pattern of decreasing call rates as the night progresses). To address these autocorrelation problems, I've tried various ARMA correlations, with correlation between time of night blocks, grouped by day (or week) and sex: gamm( .., correlation=corARMA(form=~TimeofNight | DayofYearSex), p, q) (p=1:3, q=0:3) gamm( .., correlation=corARMA(form=~TimeofNight | WeekofYearSex), p, q) This improves the ACF, but there are still correlations between residuals at same time of night - e.g. peaks every 10 lags. 1) Can I rely on the ACF when my time series isn't actually consecutive time bins, since it only includes each night (not full days). This means, for example, that the time lag between observations 9 & 10 (~1 hour) is much longer than that between 10 & 11 (~12 hours = daytime)? 2) (if yes to the above) How can I construct a correlation structure in a GAMM which allows for correlation between adjacent time bins within each night, and between the same time bin on different nights? I think you can use corCAR1, which is intended for continuous-time models, but it's limited -- it does only exponential correlations in time (analogous to corAR1). If you have enough data/computational power, you might try a 2D spline (tensor product? I don't quite remember) to allow for a gradually changing nocturnal time pattern over the course of the year (i.e something like s(timeofnight,dayofyear)) -- not entirely clear to me why you're fitting s(month) rather than s(dayofyear), unless you only have one observation per month? In general my experience is that when I find myself fitting high-order ARMA models (i.e. more than a couple of lags in total), it's because there's some systematic pattern that I failed to capture in the fixed-effects part of the model. Your mileage may of course vary. Ben Bolker ------------------------------ Message: 3 Date: Wed, 5 Jun 2013 11:45:13 +1000 From: John Maindonald <john.maindonald at anu.edu.au> To: Ben Bolker <bbolker at gmail.com> Cc: r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME] Negative Variance Message-ID: <6C1A3D6B-E7C2-4136-A511-BFEE53D983CF at anu.edu.au> Content-Type: text/plain; charset=us-ascii A variance components model that has a variance structure block variance + plot (within block) variance + subplot (within plot)
variance
makes sense only if blocks take out some part of the variation, i.e.,
variation
between plots within blocks is (in the absence of treatment effects) smaller than variation between plots in different blocks. Similarly for subplots within/between plots. If on the contrary, there is more variation between between plots within
blocks
than between plots in different blocks (this is likely to happen if there is
a
nutrient or fertility or moisture gradient within blocks), then a model that
has
the form on the second line above will if allowed account for this by
returning
a negative block component of variance estimate. It does this in order to
get
a plausible variance-covariance structure. Of course, once a gradient has been identified, it can be accommodated in the model. This does not however undo all the malign effects of an unfortunate experimental design. John Maindonald email: john.maindonald at anu.edu.au phone : +61 2 (6125)3473 fax : +61 2(6125)5549 Centre for Mathematics & Its Applications, Room 1194, John Dedman Mathematical Sciences Building (Building 27) Australian National University, Canberra ACT 0200. http://www.maths.anu.edu.au/~johnm On 05/06/2013, at 11:05 AM, Ben Bolker <bbolker at gmail.com> wrote:
John Maindonald <john.maindonald at ...> writes:
Negative variance estimates can be very useful in alerting that the variance-covariance structure is not what one expects. Or they may allow the simplest way of specifying the overall variance-covariance structure, short of specifying the variance-covariance structure in some detail. I was told of an experiment where the experimenters had chosen blocks to be at right angles to the river bank, accordingly maximising between plot variance. This came to light, in data analysed away from the scene of the original trial, when the block variance was estimated as negative -- a very useful diagnostic. Certainly, one can check on such a possibility by specifying a suitable variance-covariance structure, but how many analysts will take that trouble?
I don't quite get the geometry you're talking about, but I take the general point that diagnostics are good and that one wouldn't necessarily think to consider negative correlation.
Or one has results from each of two eyes per person. After allowing for any systematic left/right difference, are two eyes from the same individual more or less different than two eyes from different individuals? I doubt that there is a general answer that applies to all types of eye measurements.
I find this one a little bit less convincing -- here it would seem to be perfectly natural to fit a model that allowed for positive or negative correlation. The fact remains that, whether or not it's a good idea, this is very hard to do in nlme/lme4 for structural reasons. Luckily people are suggesting alternative packages. (Anyone who would like to edit http://glmm.wikidot.com/pkg-comparison accordingly is welcome to do so ...) I don't think "how do I estimate negative variances?" has quite risen to the level of a FAQ yet, so I won't bother adding it to http://glmm.wikidot.com/faq (although again, if anyone wants to take the initiative to do so I wouldn't complain). Ben Bolker
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
_______________________________________________ R-sig-mixed-models mailing list R-sig-mixed-models at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models End of R-sig-mixed-models Digest, Vol 78, Issue 11 **************************************************
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- __________________________________________________ Jonas Klasen PhD student Genome Plasticity and Computational Genetics Max Planck Institute for Plant Breeding Research