[R-sig-ME] Application of coxme

Jonas Klasen klasen at mpipz.mpg.de
Wed Jun 5 14:02:00 CEST 2013 

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

More information about the R-sig-mixed-models mailing list