n this assignment you will apply nonlinear modellingto physiological data collected on three NZ tree fern species (cyadea =Cyatheadealbata, cyasmi = Cyathea smithii, dicsqu= Dicksoniasquarrosa) . The data was collected with a fluorometerinstrument which was used to measure the electron transport rate in photosystemII (etr) (a proxy for photosynthesis) as a function of light (par).Please find the data in theprovided
etr_lc.csvdata set. There were three individuals of each species (ind) and those are true statistical replicates, so no need to aggregate the data in this case. You will see that the etrdeclines at high light levels indicating photoinhibition. Use the attached self-starter function of theEilers& Peeters(1988) photosynthesis model, whichaccounts for this photoinhibitoryeffect. Run a varying coefficient model allowing individual parameter estimates for each tree fern species. Check the model assumptions and apply variance modelling techniques if required. Then formulate a model with common parameter estimates for all species using the same variance structure as in the varying coefficient model and compare the two models using a likelihood ratio test. If the varying coefficient model is supported by the likelihood ratio test, follow up with a
post-hoccomparison using the attached custom pwcfunction (you also find the code for this function in the nonlinear regression course PDF).
pwc <- function(model,="" params,="" p.adj_method="p.adjust.methods," digits="3)" {="" #="" model=""><- gnls_vp11="" #="" params=""><- c("asym",="" "xmid")="" #="" rm(asterisks)="" mctab=""><- data.frame()="" for="" (i="" in="" seq_along(params))="" {="" p=""><- names(coef(model))[grep(pattern="params[i]," names(coef(model)))]="" p_mat=""><- contrmat(n="table(p)," type="Tukey" )="" mc=""><- apply(p_mat,="" margin="1," fun="function(x)" anova(model,="" l="x)" )="" mctab=""><- rbind(mctab,="" do.call(what="rbind," args="lapply(mc," fun="as.data.frame)))" }="" mctab=""><- cbind(mctab,="" `p-value_adj`="p.adjust(p" =="" mctab$`p-value`,="" method="ifelse(test" =="" missing(p.adj_method),="" yes="p.adj_method"><- "bh",="" no="match.arg(p.adj_method))))" asterisks=""><- ifelse(mctab$`p-value_adj`="">< 0.001,="" "***",="" ifelse(mctab$`p-value_adj`="">= 0.001 & mctab$`p-value_adj` < 0.01,="" "**",="" ifelse(mctab$`p-value_adj`="">= 0.01 & mctab$`p-value_adj` < 0.05,="" "*",="" ifelse(mctab$`p-value_adj`="">= 0.05 & mctab$`p-value_adj` <= 0.1,="" ".",="" "="" "))))="" mctab[,="" 2:ncol(mctab)]=""><- apply(mctab[,="" 2:ncol(mctab)],="" margin="2," formatc,="" digits="3," format="f" )="" mctab[as.numeric(mctab$`p-value_adj`)="">< 0.001,="" "p-value_adj"]=""><->< 0.001"="" mctab[as.numeric(mctab$`p-value`)="">< 0.001,="" "p-value"]=""><->< 0.001"="" mctab=""><- cbind(mctab,="" "="" "="format(asterisks," justify="left" ))="" mctab="" }="" ##="========================================##" ####="" eilers-peeters="" light="" response="" model="" ####="" ##="========================================##" ##="" self-starter="" funtion="" for="" the="" eilers-peeters="" light="" response="" model="" accounting="" for="" photoinhibition="" (eilers="" &="" peeters,="" 1988)="" ##="" model="" parameters:="" ##="" etrmax="maxmimum" electron="" transport="" rate="" derived="" from="" chlorophyll="" fluorescence="" ##="" measurements="" (to="" be="" replaced="" by="" psmax,="" the="" maximal="" rate="" of="" photosynthesis="" when="" ##="" dealing="" with="" gas="" exchange="" measurements)="" ##="" ppfdsat="light" saturation="" point,="" corresponds="" to="" the="" x-value="" at="" etrmax="" (ppfd="Photosynthetic" photon="" flux="" density)="" ##="" a="light-use-efficiency," initial="" slope="" of="" the="" curve="" ep=""><- function(ppfd,="" etrmax,="" ppfdsat,="" alpha)="" {="" etrmax="" *="" 2="" *="" (1="" +="" alpha)="" *="" (ppfd/ppfdsat)/((ppfd/ppfdsat)^2="" +="" 2="" *="" alpha="" *="" ppfd/ppfdsat="" +="" 1)="" }="" #="" ep=""><- function(ppfd,="" etrmax,="" ppfdsat,="" alpha)="" {="" ppfd="" (="" ppfd^2/(alpha="" *="" ppfdsat^2)="" +="" ppfd/etrmax="" -="" 2="" *="" ppfd/(alpha="" *="" ppfdsat)="" +="" (1/alpha))="" }="" ep_init=""><- function(mcall,="" lhs,="" data)="" #="" initial="" value="" routine="" for="" parameter="" estimation="" {="" xy=""><- sortedxydata(x="mCall[["ppfd"]]," y="LHS," data="data)" #="" xy=""><- sortedxydata("ppfd",="" "etr",="" data="treefern)" #="" etrmax=""><- max(xy$x)="" #="" max="" etr="" rate="" #="" ppfdsat=""><- xy$x[which.max(xy$y)]="" #="" light="" saturation="" #="" alpha=""><- coef(lm(y="" ~="" x,="" data="xy[1:3," ]))[[2]]="" #="" initial="" slope="quantum" yield="" (light="" use="" efficiency)="" sp=""><- smooth.spline(x="xy$x," y="xy$y)" #="" plot(y="" ~="" x,="" data="xy)" xv=""><- seq(min(xy$x),="" max(xy$x),="" length.out="200)" #="" lines(xv,="" predict(object="sp," x="xv)$y)" preds=""><- as.data.frame(predict(sp,="" x="xv," deriv="1))" #="" first="" derivative="" of="" the="" smoothing="" spline="" alpha=""><- mean(preds[preds$x="">< 40,="" "y"])="" #="" initial="" slope="quantum" yield="" etrmax=""><- nlsstrtasymptote(xy)="" #="" max="" etr="" rate="" ppfdsat=""><- nlsstclosestx(xy,="" yval="NLSstRtAsymptote(xy))" #="" light="" saturation="" init.pars=""><- c(etrmax,="" ppfdsat,="" alpha)="" names(init.pars)=""><- mcall[c("etrmax",="" "ppfdsat",="" "alpha")]="" init.pars="" }="" ssep=""><- selfstart(ep, ep_init, c("etrmax", "ppfdsat", "alpha")) selfstart(ep,="" ep_init,="" c("etrmax",="" "ppfdsat",="">