Temporal trend in nest survival - RMark

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

Temporal trend in nest survival - RMark

Postby cgres » Sun Apr 24, 2016 5:27 pm

Hi,

I am using RMark and have a question. I hope you can help with this.
I’ve been reading MARK and RMark books. I found in chapter 6 how to analyze a temporal trend. However, it is still not completely clear how should I write my model because all my variables are continuous, so then year is an individual covariate, not a factor.

I want to analyze if there is a change over time on the seasonal variation of nest survival and also get DSR for each year using:
- year=1,2,3,4...
- relative laying date (RLD) where the median = 0
- nest age

I tried:
A) using the original value of year as a factor for grouping (yearf=1990,1991, etc.)
B) I also substituted RLD by ~Time

See here the code for A):
Code: Select all
df.process <- process.data(nestDF,model="Nest", nocc=41, groups = "yearf")
df.ddl <- make.design.data(df.process)
full.model=list(formula=~1+RLD+year+I(year^2)+RLD:year +RLD:I(year^2+NestAge)
NS1=mark(df.process, df.ddl, model.parameters=list(S = full.model))


My design matrix doesn't show a matrix (because I have only continuous variables and the matrix is created using contrasts), and it only shows one group.
This is what I get when I try to call the matrix:
S:(Intercept) S:RLD S:year S:year2 S:NestAge S:RLD:year S:RLD:year2
S g1990 a0 t1 "1" "RLD" "year" "year2" "NestAge" "product(year,RLD)" "product(year2,RLD)"

I assume this is a major problem since there is no variation on the values for each occasion nor year.

See here the code for B):
Code: Select all
full.model=list(formula=~1+Time+yr+I(yr^2)+Time:yr +Time:I(yr^2)+NestAge)


What is the best way to write the model or how can I analyze my data?

Thanks a lot in advance,
cgres
 
Posts: 6
Joined: Wed Jan 27, 2016 6:51 pm

Re: Temporal trend in nest survival - RMark

Postby jlaake » Mon Apr 25, 2016 12:15 pm

Cynthia-

You need to provide some more details. Please contact me offlist and send me snippets of your data, design data and code you are using to do the analysis.

--jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Temporal trend in nest survival - RMark

Postby jlaake » Mon Apr 25, 2016 5:30 pm

Cynthia-

You have done a reasonable job at jumping into this but I think because you are new and haven't spent enough time with the documentation in MARK and RMark that you are missing critical pieces which is causing you confusion.

First of all, you needn't have made yr an individual covariate. In comparison reLDM is naturally an individual covariate although there are ways around that as well. First in regards to yr I changed your code to:

Code: Select all
nest = read.csv("RMarkNS16.csv", header=TRUE)
nest$yr=NULL

nocc <- max(nest$LastChecked) #number of 'daily' encounter occasions
n <- dim(nest)[1] # sample size of nests monitored
nest$Year <- as.factor(nest$year) #treat one of the year columns as a factor
nest$year=NULL

# process the data for RMark
gsgo.process <- process.data(nest,model="Nest", nocc=nocc, begin.time="1990", groups = "Year")#
gsgo.ddl <- make.design.data(gsgo.process)
gsgo.ddl$S$yr=as.numeric(as.character(gsgo.ddl$S$Year))
gsgo.ddl$S$yr=gsgo.ddl$S$yr-min(gsgo.ddl$S$yr)


Now yr is a continuous covariate in the design data. The following is a better way to set up your models and here I use the Time models.

Code: Select all
do_1=function()
{
# ~Time instead of Relative Laying Date
S.1=list(formula=~1+Time+yr+I(yr^2)+Time:yr +Time:I(yr^2)+NestAge)
S.2=list(formula=~1+Time+yr+I(yr^2) +Time:yr+NestAge)
S.3=list(formula=~Time+yr+I(yr^2)+ Time:I(yr^2)+NestAge)
S.4=list(formula=~Time+yr+I(yr^2)+NestAge)
S.5=list(formula=~Time+yr+NestAge)
S.6=list(formula=~Time+NestAge)
S.7=list(formula=~yr+NestAge)
S.8=list(formula=~NestAge)
S.9=list(formula=~yr)
S.10=list(formula=~Time)
S.11=list(formula=~Time+yr+I(yr^2)+Time:yr+Time:I(yr^2))
S.12=list(formula=~Time+yr+I(yr^2)+Time:yr)
S.13=list(formula=~Time+yr+I(yr^2))
S.14=list(formula=~Time+yr+Time:yr)
S.15=list(formula=~Time+yr)
S.16=list(formula=~1)
cml=create.model.list("Nest")
return(mark.wrapper(cml,data=gsgo.process,ddl=gsgo.ddl))
}

gsgo.results=do_1()


If I want to look at the design matrix for a model I can use:

Code: Select all
gsgo.results[[3]]$design.matrix


The row numbers in the design matrix match the PIM index value in

Code: Select all
PIMS(gsgo.results[[3]],"S")


Now if I go to your reLDM model I use the code:

Code: Select all
do_0=function()
{
# ~Time instead of Relative Laying Date
S.1=list(formula=~1+reLDM+yr+I(yr^2)+reLDM:yr +reLDM:I(yr^2)+NestAge) #full model
S.2=list(formula=~1+reLDM+yr+I(yr^2) +reLDM:yr+NestAge)
S.3=list(formula=~1+reLDM+yr+I(yr^2)+ reLDM:I(yr^2)+NestAge)
S.4=list(formula=~1+reLDM+yr+I(yr^2)+NestAge)
S.5=list(formula=~1+reLDM+yr+NestAge)
S.6=list(formula=~1+reLDM+NestAge)
S.7=list(formula=~1+yr+NestAge)
S.8=list(formula=~1+NestAge)
S.9=list(formula=~1+yr)
S.10=list(formula=~1+reLDM)
S.11=list(formula=~1+reLDM+yr+I(yr^2)+reLDM:yr+reLDM:I(yr^2))
S.12=list(formula=~1+reLDM+yr+I(yr^2)+reLDM:yr)
S.13=list(formula=~1+reLDM+yr+I(yr^2))
S.14=list(formula=~1+reLDM+yr+reLDM:yr)
S.15=list(formula=~1+reLDM+yr)
S.16=list(formula=~1)
cml=create.model.list("Nest")
return(mark.wrapper(cml,data=gsgo.process,ddl=gsgo.ddl))
}

gsgo.results_0=do_0()



And I can look at the design matrix as:

Code: Select all
> gsgo.results_0[[10]]$design.matrix
              S:(Intercept) S:reLDM S:yr S:I(yr^2) S:NestAge S:reLDM:I(yr^2)   
S g1990 a0 t1 "1"           "reLDM" "0"  "0"       "NestAge" "0"               
S g1991 a0 t1 "1"           "reLDM" "1"  "1"       "NestAge" "product(reLDM,1)"
S g1992 a0 t1 "1"           "reLDM" "2"  "4"       "NestAge" "product(reLDM,4)"
S g1993 a0 t1 "1"           "reLDM" "3"  "9"       "NestAge" "product(reLDM,9)"
S g1994 a0 t1 "1"           "reLDM" "4"  "16"      "NestAge" "product(reLDM,16)"
S g1995 a0 t1 "1"           "reLDM" "5"  "25"      "NestAge" "product(reLDM,25)"


There are now just 6 parameters:

Code: Select all
> PIMS(gsgo.results_0[[10]],"S")
group = Year1990
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
1 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
group = Year1991
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
1 2 2 2 2 2 2 2 2 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
group = Year1992
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
1 3 3 3 3 3 3 3 3 3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
group = Year1993
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
1 4 4 4 4 4 4 4 4 4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
group = Year1994
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
1 5 5 5 5 5 5 5 5 5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5
group = Year1995
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
1 6 6 6 6 6 6 6 6 6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6
>


and the estimates are:
Code: Select all
                         estimate        se        lcl        ucl
S:(Intercept)    4.5972662 0.1760142  4.2522783  4.9422541
S:reLDM         -0.1677550 0.0482490 -0.2623230 -0.0731870
S:yr             0.3553767 0.1510127  0.0593919  0.6513615
S:I(yr^2)       -0.1710913 0.0268510 -0.2237192 -0.1184634
S:NestAge        0.0284599 0.0104346  0.0080081  0.0489118
S:reLDM:I(yr^2)  0.0110755 0.0025051  0.0061656  0.0159854



You can use these with the values you want to compute the real parameters using the inverse logit link or you can use the RMark function covariate.predictions to compute the values for S for a range of reLDM values. The values for yr are one of the 6 index values.

You can use these with the individual covariate values (reLDM and NestAge) you want to compute the real parameters using the inverse logit link or you can use the RMark function covariate.predictions to compute the values for S for a range of reLDM values. The values for yr are one of the 6 index values. Index 1 is for 1990,..., and index 6 is for 1995.

I hope this helps. I didn't see anything wrong with the way you set up the models but hopefully the list of models makes sense to you. Asking whether there is a difference in the effect of reLDM and yr is simply an interaction. Currently you have the interaction parameterized as a smooth trend over time. Whether that is what you want or not depends on the situation. Alternatively, you could have used Year in place of yr and then Year*reLDM would fit 6 different lines with different intercepts and different slopes. With yr*reLDM both the intercepts and slopes are changing linearly.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Temporal trend in nest survival - RMark

Postby JClarke » Wed Jan 04, 2023 4:06 am

I am struggling with, what I suspect, is a simple question related to the original question.

When your top model has a time covariate and a second covariate that varies independently of time, how would you use covariate.predict to show how DSR varies with time within multiple groups?
I've done with in models where I have two continuous covariates and can hold one constant at its mean while I plot the other but I'm struggling with shows the impacts of time.

The model I am trying to plot is as follows:

BWTE.mod <- mark(BWTE.surv,
nocc=max(BWTE.surv$LastChecked),
model = "Nest",
groups = "cTreat",
model.parameters = list(S = list(formula = ~-1 + cTreat + Time + LitterD)))

I am relatively new to RMark and have been working my way through the documentation and workshop notes as well as the MANY phidot forums but haven't come across an answer to this yet so I would really appreciate any help that can be offered.
JClarke
 
Posts: 4
Joined: Thu Dec 29, 2022 8:37 pm

Re: Temporal trend in nest survival - RMark

Postby jlaake » Wed Jan 04, 2023 9:19 am

From your description, what you have within each group cTreat is an individual covariate LitterD which is an individual nest specific intercept and Time is a slope. Thus on the link scale this is a series of lines and something nonlinear on the real scale. What you want to do is select a few of the parameter indices across time and specify them in the indices argument as a vector. You'll find the values in the design data. Note that they will differ across groups (cTreat). Then specify your individual covariate values in data. The covariate.predictions function will compute the real value for each index with each covariate value. For any one covariate value a plot of the line is what you want and to show the effect of both you will have a series of lines. You may want to use a separate plot for each group.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Temporal trend in nest survival - RMark

Postby JClarke » Thu Jan 05, 2023 11:05 am

Thank you for clarifying this for me. For anyone who may want to do something similar, I was able to get the graphs I was looking for with the following code. It extracts indices for three different days for each treatment and plots them onto a single graph (I'll likely plot them on multiple graphs in the future.

Code: Select all
plotdata <- CCSP4.results$S.height

minheight <- min(CCSP.surv$Veg.Height)
maxheight <- max(CCSP.surv$Veg.Height)
Veg.Heightvalues = seq(from = minheight,
                       to = maxheight,
                       length = 100)

minbhco <- min(CCSP.surv$BHCONum)
maxbhco <- max(CCSP.surv$BHCONum)
BHCONum.values = seq(from = minbhco,
                     to = maxbhco,
                     length = 100)

height.pred <- covariate.predictions(plotdata,
                                     data=data.frame(Veg.Height=Veg.Heightvalues,
                                                  BHCONum=mean(BHCONum.values)),
                                     indices = c(6, 41, 81, 346, 481, 421))

Day5.21 <- which(height.pred$estimates$par.index == 6)
Day40.21 <- which(height.pred$estimates$par.index == 41)
Day80.21 <- which(height.pred$estimates$par.index == 81)
Day5.22 <- which(height.pred$estimates$par.index == 346)
Day40.22 <- which(height.pred$estimates$par.index == 481)
Day80.22 <- which(height.pred$estimates$par.index == 421)

height.pred$estimates$Date <- NA
height.pred$estimates$Date[Day5.21] <- "PBG Early"
height.pred$estimates$Date[Day40.21] <- "PBG Mid"
height.pred$estimates$Date[Day80.21] <- "PBG Late"
height.pred$estimates$Date[Day5.22] <- "SLG Early"
height.pred$estimates$Date[Day40.22] <- "SLG Mid"
height.pred$estimates$Date[Day80.22] <- "SLG Late"
head(height.pred$estimates)

library(ggplot2)

height.plot <- ggplot(height.pred$estimates,
                      aes(x = Veg.Height,
                          y = estimate,
                          group = Date)) +
  geom_line(size = 1.5,
            aes(colour = Date)) +
  geom_ribbon(aes(ymin = lcl,
                  ymax = ucl),
              alpha = 0.1) +
  xlab("Veg Height (mm)") +
  ylab("Estimated DSR") +
  ylim(0.7, 1) +
  theme_bw()
JClarke
 
Posts: 4
Joined: Thu Dec 29, 2022 8:37 pm

Re: Temporal trend in nest survival - RMark

Postby chasingbirds » Thu Jun 15, 2023 5:03 pm

Hi all,
I'm having similar issues as the folks above. My top model is an interaction effect between minimum daily temperature * elevation. I'm trying to hold elevation constant at 4 different elevations (200m, 500m, 800m and 1,200m) to see how minimum daily temperature effects DSR when elevation is held constant.

I'm trying to follow JClarke ^ above but I'm running into an error when attempting to run the covariate predictions. I'm guessing something is severly wrong with the way I'm following the code... perhaps I'm unsure what my "plotdata" should be, as well as when I'm setting the minimum and maximum values for both my covariates. Below is my code, any help is appreciated, thank you!

Code: Select all
#here is the actual model
ElevMinhypoth=mark(init.pr,ddl=init.ddl, nocc=63, model.parameters=list(S=list(formula=~elevation*min_temp)))

plotdata <- ElevMinhypoth$results

#weather3 is the dataframe with minimum daily temp in it
mintemp<- min(weather3$min_temp)
maxtemp<-max(weather3$min_temp)

tempvalues = seq(from =mintemp,
                 to=maxtemp,
                 length = 100)

#and elevation
minelev<- min(weather3$elevation)
maxelev<- max(weather3$elevation)

elev_values = seq(from=minelev,
                  to=maxelev,
                  length=100)
elev_values

#error here
temp.pred <- covariate.predictions(ElevMinhypoth2, data=data.frame(temp=tempvalues, elevation=(elev_values)), indices=c(200, 500, 800, 1200))
chasingbirds
 
Posts: 9
Joined: Wed Dec 11, 2019 1:08 pm

Re: Temporal trend in nest survival - RMark

Postby jlaake » Fri Jun 16, 2023 8:30 am

You didn't provide the error message so it is hard to discern the problem. However, you have used min_temp in the model and called the variable temp in the dataframe for covariate.predictions. The variables must have the same names and I'm guessing that is the problem.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: Temporal trend in nest survival - RMark

Postby chasingbirds » Fri Jun 16, 2023 10:58 am

Apologies... here is the code and error message below. I changed "temp" to "min_temp" and also specified certain elevations for the indices that were within "tempvalues" but still was unsuccessful. This seems like an issue with the model and the data, so I've included some info about min_temp and elevation below. Also, I tried another method to save the dataframe and then insert it into "covariate.predictions" but am still getting an error... so I'm sure this has something to do with the data.frame and the model, but I'm unsure where to edit things. In the second box of coding, I included the model results output to provide a bit more detail. Thank you!!

Code: Select all
> temp.pred <- covariate.predictions(ElevMinhypoth2, data=data.frame(min_temp=tempvalues, elevation=elev_values), indices=c(281, 507, 803, 1143))
Error in fixedvalues[fixedparms] <- model$fixed$value[match(indices[fixedparms],  :
  replacement has length zero

#here's a small example of the two covariates
> tempvalues
  [1] -2.000000000 -1.713131313 -1.426262626 -1.139393939 -0.852525253 -0.565656566
  [7] -0.278787879  0.008080808  0.294949495  0.581818182  0.868686869  1.155555556
 [13]  1.442424242  1.729292929  2.016161616  2.303030303  2.589898990  2.876767677
 [19]  3.163636364  3.450505051  3.737373737  4.024242424  4.311111111  4.597979798

> elev_values
  [1]  281.0000  289.7071  298.4141  307.1212  315.8283  324.5354  333.2424  341.9495
  [9]  350.6566  359.3636  368.0707  376.7778  385.4848  394.1919  402.8990  411.6061
 [17]  420.3131  429.0202  437.7273  446.4343  455.1414  463.8485  472.5556  481.2626


And here is another method I tried but didn't have success with
Code: Select all
##Here's the model output
> plotdata<- ElevMinhypoth2$results
> plotdata
$lnl
[1] 143.4213

$deviance
[1] 143.4213

$deviance.df
[1] 332

$npar
[1] 4

$n
[1] 535

$AICc
[1] 151.4967

$beta

$real

$beta.vcv
              ROW1          ROW2          ROW3          ROW4
[1,]  4.4557612638 -5.698525e-03 -3.618606e-01  4.537946e-04
[2,] -0.0056985249  8.494462e-06  4.528849e-04 -6.651830e-07
[3,] -0.3618605600  4.528849e-04  3.154749e-02 -3.854646e-05
[4,]  0.0004537946 -6.651830e-07 -3.854646e-05  5.538754e-08

$derived
$derived$`S Overall Survival`


$derived.vcv
$derived.vcv$`S Overall Survival`
              ROW1          ROW2         ROW3         ROW4          ROW5          ROW6
[1,]  1.001997e-04 -1.554079e-04 1.480637e-04 1.988309e-05  0.0003640541 -1.157045e-05
[2,] -1.554079e-04  6.661518e-03 1.297141e-05 1.578953e-03 -0.0033957456  2.066732e-03
[3,]  1.480637e-04  1.297141e-05 2.948564e-04 1.868153e-04  0.0009896062  3.473125e-06
[4,]  1.988309e-05  1.578953e-03 1.868153e-04 6.701941e-04  0.0015325129  4.596624e-04
[5,]  3.640541e-04 -3.395746e-03 9.896062e-04 1.532513e-03  0.0223332398 -1.083932e-03
[6,] -1.157045e-05  2.066732e-03 3.473125e-06 4.596624e-04 -0.0010839319  7.109052e-04
[7,]  1.451296e-04  3.225012e-05 3.012646e-04 2.501047e-04  0.0015609933  1.106732e-05
[8,]  1.514292e-05  1.238615e-03 1.086193e-04 5.120116e-04  0.0013362635  4.019897e-04
             ROW7         ROW8
[1,] 1.451296e-04 1.514292e-05
[2,] 3.225012e-05 1.238615e-03
[3,] 3.012646e-04 1.086193e-04
[4,] 2.501047e-04 5.120116e-04
[5,] 1.560993e-03 1.336264e-03
[6,] 1.106732e-05 4.019897e-04
[7,] 3.261340e-04 1.674928e-04
[8,] 1.674928e-04 4.247997e-04


$covariate.values
NULL

$singular
NULL

$real.vcv
NULL

#creating the dataframe
predict1<-data.frame(min_temp=seq(min(weather3$min_temp), max(weather3$min_temp), length.out=100), elevation=seq(min(weather3$elevation), max(weather3$elevation), length.out=50))

min_temp          elevation
-2.000000000   281.0000         
-1.713131313   298.5918         
-1.426262626   316.1837         
-1.139393939   333.7755   


#writing the covariate.predictions model
> temp.pred2<-covariate.predictions(ElevMinhypoth2, data=predict1)$estimates
Error in deriv.real %*% model$results$beta.vcv :
  non-conformable arguments


chasingbirds
 
Posts: 9
Joined: Wed Dec 11, 2019 1:08 pm

Re: Temporal trend in nest survival - RMark

Postby jlaake » Mon Jun 19, 2023 12:27 pm

The problem here was that there were only 496 values of model.index and some values specified were greater than 496 because of a misunderstanding of what indices specified. I need to add some code to catch this error and provide a more informative error message.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to RMark

Who is online

Users browsing this forum: No registered users and 14 guests

cron