[Basta-users] Calculating daily mortality rates for models with cont. covariates using estimates
Fernando Colchero
colchero at imada.sdu.dk
Tue Mar 29 10:22:57 CEST 2016
Hi Dries,
Definitely, that will be correct, add the other parameters and take a fixed value for the covariates that you are not showing. I agree that taking the mean would make more sense and it would make it clearer to explain. Your updated mortality is correct.
Best,
Fernando
Fernando Colchero
Assistant Professor
Department of Mathematics and Computer Science
Max-Planck Odense Center on the Biodemography of Aging
Tlf. +45 65 50 23 24
Email colchero at imada.sdu.dk <mailto:colchero at imada.sdu.dk>
Web www.sdu.dk/staff/colchero <http://www.sdu.dk/staff/colchero>
Pers. web www.colchero.com <http://www.sdu.dk/staff/colchero>
Adr. Campusvej 55, 5230, Odense, Dk
University of Southern Denmark
> On 28 Mar 2016, at 15:22, Dries Van de Loock <dries.vandeloock at ugent.be> wrote:
>
> Dear Fernando,
>
> many, many thanks for this. It was exactly what I was looking for and will greatly improve our manuscript.
>
> I have one questions remaining.
>
> How should I approach this method when wanting to visualise the effect of one continuous variable from a multivariate model both including continuous and one categorical variable with 5 levels ?
>
> - Could I use the weighted average (=average c, b1 and b0 based on cohorts sample sizes) of the mortality estimates ? When just averaging, I have found that the mort rates were far too high as one cohort with just one case had high estimates.
> - Could I just add another term * exp(g * z) to the model for each continuous variable and use the averaged predictor (z) value for the variable I'm not interested in ?
>
>
> For example, we have a model : mort ~ LOC + SMI + GS
>
> GS = groupsize (continuous)
> SMI = mass (continuous, like above example)
> LOC = location (5 levels, sample sizes differ hugely between cohorts)
>
> When I want to visualise the effect of SMI solely, would this approach be correct :
>
> #define function
> mort <- function(x, z1, z2, b, g) {
> exp(g[1] * z1) * exp(g[2] * z2) * (b[1] + b[2] * b[3]^(b[2]) * x^(b[2] - 1))
> }
>
> # define values
> x <- seq(1, 55, 1)
> z1 <- c(min(covariates$SMI),max(covariates$SMI)) # min and max value of the continuous variable of interest
> z2 <- mean(covariates$GS) # mean of the continuous variable not interested in to visualise
> c <- mean(out$coefficients[c(1:5),1]) # average of the c values of each LOC (c.LOC1, c.LOC2, etc.) (in this example not yet weighted)
> b0 <- mean(out$coefficients[c(6:10),1]) # average of the b0 values of each LOC (b0.LOC1, b0.LOC2, etc.) (in this example not yet weighted)
> b1 <- mean(out$coefficients[c(11:15),1])
> b <- c(c,b0,b1)
> g <- out.SMI$coefficients[c(16:17), 1] # g parameters
>
> # calculate
> mumin <- mort(x, z1[1], z2, b, g)
> mumax <- mort(x, z1[2], z2, b, g)
>
>
> I hope the example and my question are clear.
>
> Again, many thanks for your time invested in this and I'm looking forward to your reply !
>
> Dries
>
>
> 2016-03-26 8:25 GMT+01:00 Fernando Colchero <colchero at imada.sdu.dk <mailto:colchero at imada.sdu.dk>>:
> Dear Dries,
>
> Sorry for the late reply. This is certainly an issue we have not yet addressed with the package. Alternatively, you could produce a plot by extracting the information from the coefficients table in the BaSTA output. Here’s an option based on your results:
>
> # Define the Weibull mortality with proportional hazards as:
> mort <- function(x, z, b, g) {
> exp(g * z) * (b[1] + b[2] * b[3]^(b[2]) * x^(b[2] - 1))
> }
>
> where x is age, z is the SMI value, b is the vector of mortality parameters and g is the vector of proportional hazards parameters.
>
> # Define x (I’m making up this range of ages…):
> x <- seq(0.1, 5, 0.1)
>
> # Define the lower and upper bounds of z (these are also made up…):
> z <- c(3, 6)
>
> # Extract the b and g parameters:
> b <- out$coefficients[1:3, 1]
> g <- out$coefficients[4, 1]
>
> # Calculate the corresponding mortalities:
> mulow <- mort(x, z[1], b, g)
> muup <- mort(x, z[2], b, g)
>
> # Define the range for the y-axis:
> ylim <- c(0, max(c(mulow, muup)))
>
> # Plot the results:
> plot(x, mulow, col = 1, type = 'l', ylim = ylim)
> lines(x, muup, col = 2)
> legend('topright', c("Low SMI", "High SMI"), col = c(1, 2),
> lwd = 2)
>
> Let me know if this is of any help.
>
> Best,
>
> Fernando
>
> Fernando Colchero
> Assistant Professor
> Department of Mathematics and Computer Science
> Max-Planck Odense Center on the Biodemography of Aging
>
> Tlf. +45 65 50 23 24 <tel:%2B45%2065%2050%2023%2024>
> Email colchero at imada.sdu.dk <mailto:colchero at imada.sdu.dk>
> Web www.sdu.dk/staff/colchero <http://www.sdu.dk/staff/colchero>
> Pers. web www.colchero.com <http://www.sdu.dk/staff/colchero>
> Adr. Campusvej 55, 5230, Odense, Dk
>
> University of Southern Denmark
>
>> On 22 Mar 2016, at 11:32, Dries Van de Loock <dries.vandeloock at ugent.be <mailto:dries.vandeloock at ugent.be>> wrote:
>>
>> Dear BaSTA community,
>>
>> Thanks to your help and BaSTA, we could statistically confirm some of our expectations on post-fledging survival in a afrotropical passerine <http://lists.r-forge.r-project.org/pipermail/basta-users/2015-October/000131.html>.
>>
>> While plotfancybasta() and plot() produce nice graphs for categorical variables, we are looking for a way to visualise the impact of continues covariates on survival of the fledglings.
>>
>> More specifically, we fitted a weibull survival curve with makeham shape + mass at fledging (SMI) as a covariate (base model + 1 covariate) and find a negative relation with mortality (CI does not include zero).
>>
>> Coefficients:
>> Estimate StdErr Lower95%CI Upper95%CI SerAutocor UpdateRate PotScaleReduc
>> c 0.02704 0.01595 0.009527 0.066033 0.56051 0.2223 1.017
>> b0 0.06541 0.06751 0.001576 0.252938 0.38508 0.2605 1.002
>> b1 0.90422 0.63824 0.043231 2.358846 0.06547 0.2501 1.001
>> gamma.SMI -0.03971 0.01978 -0.083632 -0.006801 0.55284 0.2526 1.003
>> pi.1 0.14211 0.01089 0.121642 0.163994 0.03578 1.0000 1.001
>>
>> We see two ways to visualise this
>>
>> 1. Daily mortality rates over a two month period (y-axis) plotted against time (x-axis) for a discrete number of SMI values. As such, you visualise the mortality shape for a number of SMI values over a two month period.
>>
>> 2. The cumulative mortality probabilities (y-axis) after a two month period plotted against SMI (x-axis). We then expect that cumulative mortality probabilities will decrease with increasing SMI.
>>
>> To do this, we could use the estimates to calculate (i) the daily mortality rates over a two month period and use these to calculate (ii) the cumulative mortality probabilities after a two month period for a discrete number of SMI values within the measured range (for example : 20, 25,30,35).
>>
>>
>> While we successfully calculate daily mortality rates fitting the estimates in the function (b0 * b1 * ((b1 * k)^(b0 - 1))) + c) when no covariates were added to the model, we failed to find a method on how to include the gamma estimate in the function.
>>
>>
>> Any help on how we can calculate daily mortality rates for a specific value of the covariate and be able to produce these graphs is highly appreciated.
>>
>> Best wishes,
>>
>> Dries Van de Loock
>> PhD student
>> Terrestrial Ecology Unit
>> Department of Biology
>> Ghent University
>> KL Ledeganckstraat 35
>> B-9000 Ghent, Belgium
>> Phone: +32 (0)9 265 50 39 <tel:%2B32%20%280%299%20265%2050%2039>
>> http://www.ecology.ugent.be/terec/ <http://www.ecology.ugent.be/terec/>
>>
>>
>>
>> _______________________________________________
>> Basta-users mailing list
>> Basta-users at lists.r-forge.r-project.org <mailto:Basta-users at lists.r-forge.r-project.org>
>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/basta-users <https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/basta-users>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/basta-users/attachments/20160329/a3d726c7/attachment-0001.html>
More information about the Basta-users
mailing list