Mortality-growth Analyses functions -

# Mortality-growth Analyses functions - ## partial for Yi-Lin ## 宜霖的筆記 --- https://rdrr.io/github/forestgeo/ctfs/man/ --- --- ## 5.1 Mortality at Slow Growth Function and errorbars 1. function to estimate annualized mortality for the slowest growing 25% of stems in each species for data with 3 censuses 2. growth is the growth rate in the first 2 censuses 3. mort is the 0,1 variable for mortality in the 3rd census 4. needs CTFS R package --- :::success slow.grow.mort.onesp ::: :::info mortality.calculation(N, S, meantime) N : The number of individuals alive at the outset. S : Number of survivors meantime : time interval ::: --- ```r= slow.grow.mort.onesp = function(dat, quant=0.25) { slow.gr = subset(dat, dat$growth <= quantile(dat$growth, probs=quant, na.rm=T) ) Mort.Rate = mortality.calculation(N=dim(slow.gr)[1], S=dim(slow.gr[slow.gr$mort==0,])[1], meantime=mean(slow.gr$time.mort)) return(Mort.Rate) } ``` --- 1. Fn to draw an axis with transformed values, either raised to an exponent, or logged (if powerexp is NULL). 2. value is on UNtransformed scale, but it will be plotted with the (untransformed) value label at the transformed location on the axis. :::spoiler 製圖程序(右圖，error bar plot) ```r=+ ! # from Hmisc #errorbar errbar = function (x, y, yplus, yminus, cap = 0.015, main = NULL, sub = NULL, xlab = as.character(substitute(x)), ylab = if (is.factor(x) || is.character(x)) "" else as.character(substitute(y)), add = FALSE, lty = 1, type = "p", ylim = NULL, lwd = 1, pch = 16, errbar.col = par("fg"), Type = rep(1, length(y)), ...) { if (is.null(ylim)) ylim <- range(y[Type == 1], yplus[Type == 1], yminus[Type == 1], na.rm = TRUE) if (is.factor(x) || is.character(x)) { x <- as.character(x) n <- length(x) t1 <- Type == 1 t2 <- Type == 2 n1 <- sum(t1) n2 <- sum(t2) omai <- par("mai") mai <- omai mai[2] <- max(strwidth(x, "inches")) + 0.25 par(mai = mai) on.exit(par(mai = omai)) plot(NA, NA, xlab = ylab, ylab = "", xlim = ylim, ylim = c(1, n + 1), axes = FALSE, main = main, sub = sub, ...) axis(1) w <- if (any(t2)) n1 + (1:n2) + 1 else numeric(0) axis(2, at = c(seq.int(length.out = n1), w), labels = c(x[t1], x[t2]), las = 1, adj = 1) points(y[t1], seq.int(length.out = n1), pch = pch, type = type, ...) segments(yplus[t1], seq.int(length.out = n1), yminus[t1], seq.int(length.out = n1), lwd = lwd, lty = lty, col = errbar.col) if (any(Type == 2)) { abline(h = n1 + 1, lty = 2, ...) offset <- mean(y[t1]) - mean(y[t2]) if (min(yminus[t2]) < 0 & max(yplus[t2]) > 0) lines(c(0, 0) + offset, c(n1 + 1, par("usr")[4]), lty = 2, ...) points(y[t2] + offset, w, pch = pch, type = type, ...) segments(yminus[t2] + offset, w, yplus[t2] + offset, w, lwd = lwd, lty = lty, col = errbar.col) at <- pretty(range(y[t2], yplus[t2], yminus[t2])) axis(side = 3, at = at + offset, labels = format(round(at, 6))) } return(invisible()) } if (add) points(x, y, pch = pch, type = type, ...) else plot(x, y, ylim = ylim, xlab = xlab, ylab = ylab, pch = pch, type = type, ...) xcoord <- par()$usr[1:2] smidge <- cap * (xcoord[2] - xcoord[1])/2 segments(x, yminus, x, yplus, lty = lty, lwd = lwd, col = errbar.col) if (par()$xlog) { xstart <- x * 10^(-smidge) xend <- x * 10^(smidge) } else { xstart <- x - smidge xend <- x + smidge } segments(xstart, yminus, xend, yminus, lwd = lwd, lty = lty, col = errbar.col) segments(xstart, yplus, xend, yplus, lwd = lwd, lty = lty, col = errbar.col) return(invisible()) } add.alpha <- function(col, alpha=1){ if(missing(col)) stop("Please provide a vector of colours.") apply(sapply(col, col2rgb)/255, 2, function(x) rgb(x[1], x[2], x[3], alpha=alpha)) } ``` ::: --- power-transform pospower is to deal with raising negative numbers to a power fractional of 1 :::success transformAxes ::: ```r=+ ! pospower = function (x, expon) { result = sign(x) * (abs(x)^(expon))#sign回傳該值的正or負號 return(result) } transformAxes=function(value,powerexp=0.45,whichaxis=2,multiplier=6, CA=0.5) { label=as.character(value) if(!is.null(powerexp)) transformed=pospower(value,powerexp) else transformed=log(value)+multiplier if(!is.null(whichaxis)) axis(whichaxis,at=transformed,label=label, cex.axis=CA) return(list(pos=transformed,lab=label)) } ``` --- ## 5.2 Make mortality-prior growth data sets function to take data sets from 3 censuses and merge them into one file for analyses * dbh is converted to cm if the unit is mm * species are limited to those with ID level == species in the species table * gr.transform can be "log" or "power" * if log, then growth rates <= zero are eliminated * if power, then 0.45 power is used, and zero and negative rates are retained * No capital letters in plot.name :::success make.mort.gr.data ::: ```r= ! make.mort.gr.data = function(census1=fushan.tree1, census2=fushan.tree2, census3=fushan.tree3, dbhunit="cm", gr.transform="power", plot.name="fushan", species.list=fushan.spptable) { dat123a=data.frame(tag=census1$tag, sp=census1$sp, gx=census1$gx, gy=census1$gy, dbh1=census1$dbh, dbh2=census2$dbh, dbh3=census3$dbh, dd1=census1$date/365.25+1960, dd2=census2$date/365.25+1960, dd3=census3$date/365.25+1960, pom1=census1$pom, pom2=census2$pom, pom3=census3$pom, DFstatus1=census1$DFstatus, DFstatus2=census2$DFstatus, DFstatus3=census3$DFstatus) dat123a$pom1=as.numeric(as.character(dat123a$pom1)) dat123a$pom2=as.numeric(as.character(dat123a$pom2)) dat123a$pom3=as.numeric(as.character(dat123a$pom3)) # convert dbh to cm if the data are in mm if (dbhunit=="mm") { dat123a$dbh1=dat123a$dbh1/10 dat123a$dbh2=dat123a$dbh2/10 dat123a$dbh3=dat123a$dbh3/10 } # limit dbh to appropriate minimum value (1 cm) dat123b = subset(dat123a, dbh1 >= 1 & dbh2 >= 1) print(range(dat123b$dbh1, na.rm=T)) # subset species list to good species (vector of names only), which are those with IDlevel=="species" if(any(colnames(species.list)=="IDlevel")) good.spp = subset(species.list, IDlevel=="species")$sp dat123 = dat123b[dat123b$sp %in% good.spp,] print(dim(dat123)) print(head(dat123)) # Calculations for Growth in census 1-2 and death in census 2-3 # subset data to get cohort to be used for calculating growth in census 1-2 and exclude NA's live=dat123$DFstatus1=="alive" & dat123$DFstatus2=="alive" & dat123$dbh1 >= 1 gdata=subset(dat123, live) print(dim(gdata)) # find the records that are good for calculating growth gooddate = gdata$dd1>0 & gdata$dd2>0 samepom = gdata$pom1 == gdata$pom2 hasdbh = !is.na(gdata$dbh1) & !is.na(gdata$dbh2) good2 = gooddate & samepom & hasdbh # subset gdata by above conditions gdata<-subset(gdata, good2) print(dim(gdata)) # calculate growth parameters and unite in dataframe # note that dbh is in cm and time in the data is in fractional years time <-(gdata$dd2-gdata$dd1) growth <-(gdata$dbh2-gdata$dbh1)/time pctred <-gdata$dbh2/gdata$dbh1 gdata <-data.frame(gdata, time=time, growth=growth, pctred=pctred) # eliminate bad growth trees (outliers in both negative and postive direction) # these match the default settings in trimgrowth in CTFSRpackage except poms must match exactly maxgrow = 7.5 intercept = 0.9036 slope = 0.006214 err.limit = 4 if (dbhunit=="cm") intercept = intercept/10 stdev.dbh1 = slope*gdata$dbh1 + intercept good.neggrow = gdata$dbh2 > (gdata$dbh1-err.limit*stdev.dbh1) good.posgrow = gdata$growth < maxgrow good = good.neggrow & good.posgrow mortgrow=subset(gdata, good) # make mort in census 3 into 0/1 variable (I think glm() needs this) # only use alive and dead, ignoring missing and lost_stem status codes mortgrow$mort = rep(NA, dim(mortgrow)[1]) mortgrow$mort[mortgrow$DFstatus3 == "alive"]=0 mortgrow$mort[mortgrow$DFstatus3 == "dead"]=1 mortgrow$time.mort <-(mortgrow$dd3-mortgrow$dd2) # LOG transformation of growth - subset data for positive growth and make logged variables if (gr.transform == "log") { mortgrow_123 = subset(mortgrow, growth>0 & !is.na(mort)) mortgrow_123$trans_growth = log(mortgrow_123$growth) mortgrow_123$log_dbh1 = log(mortgrow_123$dbh1) mortgrow_123$log_dbh2 = log(mortgrow_123$dbh2) } # apply 0.45 power transformation to growth and retain negative and zero growth rates # log transform dbh if (gr.transform == "power") { mortgrow_123 = subset(mortgrow, !is.na(mort)) mortgrow_123$trans_growth = rep(NA, length(mortgrow_123$growth)) mortgrow_123$trans_growth[mortgrow_123$growth >= 0] = mortgrow_123$growth[mortgrow_123$growth >= 0]^0.45 mortgrow_123$trans_growth[mortgrow_123$growth < 0] = -(abs(mortgrow_123$growth[mortgrow_123$growth < 0])^0.45) mortgrow_123$log_dbh1 = log(mortgrow_123$dbh1) mortgrow_123$log_dbh2 = log(mortgrow_123$dbh2) } # add dbh category dcat = cut(mortgrow_123$dbh1, breaks=c(1,5,10,20, 10000), include.lowest=T, right=F) mortgrow_123$dcat = dcat # subset to just stems with 0 or 1 for mort (no NA) mortgrow_123_final = subset(mortgrow_123, !is.na(mort)) return(mortgrow_123_final) } ``` ### Example Usage fushan_040914=make.mort.gr.data(census1=fushan.tree1, census2=fushan.tree2, census3=fushan.tree3, dbhunit="cm", gr.transform="power", plot.name="fushan", species.list=fushan.spptable) --- ## 5.3 Fit within-species mortality-growth-dbh glms * function to fit three glms of mortality in census 2-3 as *a function of growth in the census 1-2 and size at start of census 2 for a minimum species abundance of min_N stems * if you want to run the function for different size classes, then subset the data first * the function filters out species not meeting the minimum sample size criteria * note that several different models were fit in very early exploratory analyses for the project, and were retained here, * but the final model used was UU glms are: 1. model UU = uncentered, unstandardized trans growth and logdbh are NOT standardized - so intercept is mortality at trans growth (and growth) = 0 and log dbh = 0 (dbh = 1 cm) 2. model US = BOth uncentered, Both standardized/scaled trans growth and logdbh are scaled by dividing by the standard deviation of each variable for the species, so intercept is mortality at trans growth (and growth) = 0 and log dbh = 0 (dbh = 1 cm) (same as res_glm) * the difference is that the variables are scaled by the sd, so the slopes will be on the same sd scale across species, which is not desirable * the slope will be the change in mortality probability for a one unit SD unit change in trans growth or logdbh 3. model C20S = centered on 20th quantile of trans growth, standardized/scaled log dbh trans growth centered on 20% quantile of trans growth (the first non-zero quantile) & scaled by dividing by the standard deviation of logdbh for species, so intercept is mortality at low growth for species and mean logdbh for the species 4. model C20U = centered on 20th quantile of trans growth, UNstandardized/scaled log dbh trans growth centered on 20% quantile of trans growth & UNstandardized logdbh so intercept is mortality at low growth for species and 1 cm (logdbh=0--> dbh=1 cm) requires data output from make.mort.gr.data and function and CTFS Rpackage must be loaded transformed growth rate can be either log or power, depending on what was chosen in make.mort.gr.data() also produces a PDF of the graphs of the individual fits for each species, with mean mortality rate printed in corner (on left panel) and mean mortality in growth classes on right panel * y axis is back-transformed to the probability scale (0 to 1); dbh is set at mean log dbh; points are scaled to dbh * the output of the function is a list with 4 elements * coef_glm is a data frame of the glm coefficients for all 3 models in columns for each species (in rows) * res_glm_list is a list of the model UU glm objects for all species, with the $data element of the list set to NULL (for memory size reasons); this can be used for plotting and prediction 1. res_glm_std_list is same as res_glm_list , but for model US 2. res_glm_low3_list is same as res_glm_list , but for model C20S 3. res_glm_low4_list is same as res_glm_list , but for model C20U this takes a few minutes to run, depending on how many species there are :::success mort.grow.glm.plot ::: ```r= ! mort.grow.glm.plot = function (data.to.use=lambir_970308, min_N=200, min_N_mort=5) { res_glm_low3_list = list() res_glm_low4_list = list() res_glm_std_list = list() res_glm_list = list() coef_glm = NULL coef_glm_num = NULL coef_glm_char = NULL splist2=names(which(table(data.to.use$sp)>=min_N)) data.to.use_minN1 = data.to.use[data.to.use$sp %in% splist2,] sab=as.data.frame.matrix(table(as.character(data.to.use_minN1$sp), data.to.use_minN1$mort)) splist=rownames(sab[which(sab[,2] >= min_N_mort),]) data.to.use_minN = data.to.use_minN1[data.to.use_minN1$sp %in% splist,] xlim.range = c(min(data.to.use_minN$trans_growth), max(data.to.use_minN$trans_growth)) ``` ```r= ! pdf("Mort Growth Plots by Sp.pdf", width = 11, height = 5, onefile=T, pointsize=12) for (i in 1:length(splist) ) { print(paste(i, "of ", length(splist))) spdat = subset(data.to.use_minN, sp==splist[i]) spdat$std_logd = spdat$log_dbh2/sd(spdat$log_dbh2) spdat$std_transgr = spdat$trans_growth/sd(spdat$trans_growth) # spdat$low_transgr = (spdat$trans_growth - quantile(spdat$trans_growth, prob=0.05)) # this is essentially unstandardized for most species bc 5th quantile = 0 spdat$low_transgr = (spdat$trans_growth - quantile(spdat$trans_growth, prob=0.2)) xlim.range.lowgr = c(min(spdat$low_transgr), max(spdat$low_transgr)) # model UU = uncentered, unstandardized # trans growth and logdbh are NOT standardized - so intercept is mortality at trans growth (and growth) = 0 and log dbh = 0 (dbh = 1 cm) # original model fit in MS res_glm = glm(mort ~ log_dbh2 + trans_growth, family = "binomial", data = spdat, na.action = na.omit, epsilon = 1e-8, maxit = 1000) # model US = BOth uncentered, Both standardized/scaled # trans growth and logdbh are scaled by dividing by the standard deviation of each variable for the species - # so intercept is mortality at trans growth (and growth) = 0 and log dbh = 0 (dbh = 1 cm) (same as res_glm) # the difference is that the variables are scaled by the sd, so the slopes will be comparable across species # the slope will be the change in mortality probability for a one unit SD unit change in trans growth or logdbh res_glm_std = glm(mort ~ std_logd + std_transgr, family = "binomial", data = spdat, na.action = na.omit, epsilon = 1e-8, maxit = 1000) # model C20S = centered on 20th quantile of trans growth, standardized/scaled log dbh # trans growth centered on 20% quantile of trans growth (the first non-zero quantile) & scaled by dividing by the standard deviation of logdbh for species - # so intercept is mortality at low growth for species and mean logdbh for the species res_glm_lowgr3 = glm(mort ~ std_logd + low_transgr, family = "binomial", data = spdat, na.action = na.omit, epsilon = 1e-8, maxit = 1000) # model C20U = centered on 20th quantile of trans growth, UNstandardized/scaled log dbh # trans growth centered on 20% quantile of trans growth & UNstandardized logdbh - # so intercept is mortality at low growth for species and 1 cm (logdbh=0--> dbh=1 cm) res_glm_lowgr4 = glm(mort ~ log_dbh2 + low_transgr, family = "binomial", data = spdat, na.action = na.omit, epsilon = 1e-8, maxit = 1000) res_glm_list[[i]] = res_glm res_glm_std_list[[i]] = res_glm_std res_glm_low3_list[[i]] = res_glm_lowgr3 res_glm_low4_list[[i]] = res_glm_lowgr4 res_glm_list[[i]]$data = NULL res_glm_std_list[[i]]$data = NULL res_glm_low3_list[[i]]$data = NULL res_glm_low4_list[[i]]$data = NULL # plots par(mfcol=c(1,2)) # fitted relationship: uncentered data; this is useful for visualization (centered version just slided the curve along x-axis) plot(spdat$trans_growth, fitted.values(res_glm), cex=spdat$dbh2^(1/4), main=paste(splist[i], ": Uncentered"), ylim=c(0,1), xlim=xlim.range, ylab = "Intercensus Mortality Probability", xlab = "Transformed Growth Rate (cm/y)") text(0.7*xlim.range[2], 0.8, paste(round(mean(fitted.values(res_glm)),4)) ) x=seq(min(spdat$trans_growth), max(spdat$trans_growth), by=0.02) points(x, predict(res_glm, newdata = data.frame(log_dbh2 = rep(mean(spdat$log_dbh2),length(x)), trans_growth=x), type = "response"), type="l", col="red", lwd=2) rug(jitter(spdat$trans_growth[spdat$mort==1]), side=3, col="#60606060") rug(jitter(spdat$trans_growth[spdat$mort==0]), side=1, col="#60606060") abline(h=0,col="gray20") abline(h=1,col="gray20") ##### not part of plotting #### # calculate mortality values to include in data summary table # note: the 5-y mortality probability (mean_mort) is the proportion of stems dying, # which is the maximum likelihood estimator for the binomial probabilty # the annual mortality probability (Amean_mort) is Rick's estimator Ntot = sum(table(spdat$mort) ) Stot = table(spdat$mort)[1] Dtot = table(spdat$mort)[2] mean_mort = Dtot/Ntot Amean_mort = (log(Ntot)-log(Stot))/mean(spdat$dd3 - spdat$dd2) mort.slow.gr = slow.grow.mort.onesp(dat = spdat, quant=0.25) std.mort = predict(res_glm, newdata = data.frame(log_dbh2 = 1.6, trans_growth = -3) ) coef_glm_num = rbind(coef_glm_num, c(as.numeric(res_glm$coefficients), as.numeric(res_glm_std$coefficients), as.numeric(res_glm_lowgr3$coefficients), as.numeric(res_glm_lowgr4$coefficients), mean(spdat$growth), quantile(spdat$growth, probs=0.95, na.rm=T), mean_mort, Amean_mort, mort.slow.gr$rate, std.mort, Ntot, Stot, mean(spdat$trans_growth), sd(spdat$trans_growth), mean(spdat$log_dbh2), sd(spdat$log_dbh2) ) ) coef_glm_char = rbind(coef_glm_char, as.character(spdat$sp[1]) ) ###### # back to plotting # calculate mortality rate in growth rate classes for plotting lgr_quant = quantile(spdat$trans_growth, c(0, 0.25,0.5,0.75, 1)) # note for some species, there are so many zero growth rates that 0.25 and 0.5 quantiles are both zero if ( any(table(lgr_quant)>1) ) lgr_quant = quantile(spdat$trans_growth, c(0, 0.5,0.75, 1)) lgr_class = cut(spdat$trans_growth, lgr_quant, include.lowest=T, right=F) spdat$lgr_class = lgr_class if ( length(table(spdat$mort))==1 ) { plot(1,1, main="All dead or alive") } else { # mort rates in growth categories only used for plotting "raw" data in PDF to check fits S= table(spdat$lgr_class, spdat$mort)[,1] N= table(spdat$lgr_class, spdat$mort)[,2] + table(spdat$lgr_class, spdat$mort)[,1] time.mort = mean(spdat$time.mort) lower.ci=find.climits(N,(N-S),kind="lower") upper.ci=find.climits(N,(N-S),kind="upper") mort.rate=(log(N)-log(S)) upper.rate=(log(N)-log(N-upper.ci)) lower.rate=(log(N)-log(N-lower.ci)) mort.rate[S==0]=upper.rate[S==0]=Inf upper.rate[upper.ci==N]=Inf lower.rate[lower.ci==N]=0 mort.rate[N==0]=lower.rate[N==0]=upper.rate[N==0]=NA errbar(x=as.numeric(lgr_quant[1:length(lgr_quant)-1]), y = as.numeric(mort.rate), yplus = as.numeric(upper.rate), yminus = as.numeric(lower.rate), xlab = "Beginning of log growth quartile", ylab="5-y Mortality prob.", lwd=2, pch=19, cex=1.2, main = splist[i], ylim=c(0,1), xlim=xlim.range) text(0.7*xlim.range[2], 0.8, paste(round(mean_mort,4)) ) } } dev.off() #### plotting done names(res_glm_list)=splist names(res_glm_std_list)=splist names(res_glm_low3_list)=splist names(res_glm_low4_list)=splist coef_glm = data.frame(coef_glm_num, coef_glm_char) colnames(coef_glm) = c("intercept","slope_log_dbh","slope_trans_growth", "std_intercept","std_slope_log_dbh","std_slope_trans_growth", "lowgr3_intercept","lowgr3_slope_log_dbh","lowgr3_slope_trans_growth", "lowgr4_intercept","lowgr4_slope_log_dbh","lowgr4_slope_trans_growth", "mean_growth", "q95_growth", "y5_mean_mort", "annual_mort", "mort_slow_gr", "std_mort", "N", "S", "mean_trans_growth", "sd_trans_growth", "mean_log_dbh2", "sd_log_dbh2", "sp") coef_glm$slope_trans_growth = as.numeric(coef_glm$slope_trans_growth) coef_glm$slope_log_dbh = as.numeric(coef_glm$slope_log_dbh) coef_glm$intercept = as.numeric(coef_glm$intercept) coef_glm$mean_growth = as.numeric(coef_glm$mean_growth) coef_glm$y5_mean_mort = as.numeric(coef_glm$y5_mean_mort) coef_glm$N = as.numeric(coef_glm$N) coef_glm$std_mort = as.numeric(coef_glm$std_mort) coef_glm$std_slope_trans_growth = as.numeric(coef_glm$std_slope_trans_growth) coef_glm$std_slope_log_dbh = as.numeric(coef_glm$std_slope_log_dbh) coef_glm$std_intercept = as.numeric(coef_glm$std_intercept) coef_glm$lowgr3_slope_trans_growth = as.numeric(coef_glm$lowgr3_slope_trans_growth) coef_glm$lowgr3_slope_log_dbh = as.numeric(coef_glm$lowgr3_slope_log_dbh) coef_glm$lowgr3_intercept = as.numeric(coef_glm$lowgr3_intercept) coef_glm$lowgr4_slope_trans_growth = as.numeric(coef_glm$lowgr4_slope_trans_growth) coef_glm$lowgr4_slope_log_dbh = as.numeric(coef_glm$lowgr4_slope_log_dbh) coef_glm$lowgr4_intercept = as.numeric(coef_glm$lowgr4_intercept) return(list(coef_glm=coef_glm, res_glm_list=res_glm_list, res_glm_std_list=res_glm_std_list, res_glm_low3_list=res_glm_low3_list, res_glm_low4_list=res_glm_low4_list) ) } # Example Usage fushan_glm_result_040914=mort.grow.glm.plot(data.to.use=fushan_040914, min_N=200) fushan_glm_coef_040914 = fushan_glm_result_040914$coef_glm ######## # between species relationship cor.test(log(fushan_glm_coef_040914$q95_growth), log(inv.logit(fushan_glm_coef_040914$intercept))) ``` --- data問題： 1. 沒有branch資料?!->檢查fushan.tree1 ```r= fushan.tree1$tag[duplicated(fushan.tree1$tag)] #character(0) ``` 如何加入 lai?可以要直接使用這個code? -->先檢查老師code是否都合理處理 pdf. 每種一個值代表甚麼?mean mort 左圖的上下bar代表? rug plots位於哪裡? 右邊的error bar 圖要如何解讀? x軸Beginning of log growth quartile什麼意思?四個點為四個quantile 解開function進行運算? 需要抓出glm coef--> 每種equation參數畫在圖上好抓出哪些是用CTFS package 可以用自己寫的glm loop 了解pseudoR、AIC 檢查missing stem問題 5.3 C20S(std_logd)、C20U(log_dbh2)，差在哪? 為何是20th quantile? 為何最後用UU?還是C20U? US、C20U、C20S後續分析有何用 std.mort = predict(res_glm, newdata = data.frame(log_dbh2 = 1.6, trans_growth = -3) ) 為何是1.6,-3?? 找出1.及3.的兩個處理處理1.-->5.3 min_N=200,min_N_mort=5 ![](https://i.imgur.com/E43ogcR.jpg) ```graphviz digraph summary{ f[label="slow.grow.mort.onesp"] pospower[shape=box] transformAxes[shape=box] pospower->transformAxes } ``` ```graphviz digraph summary{ start [label="fushan.tree1-3"] dat123a spdat [label="spdat"] UU [label="res_glm (UU)",shape=box] US [label="res_glm_std (US)",shape=box] C20S [label="res_glm_lowgr3 (C20S)",shape=box] C20U [label="res_glm_lowgr4 (C20U)",shape=box] mor1[label="mort.slow.gr"] mor2 [label="slow.grow.mort.onesp", shape=box] fu[label="fushan_040914"] fu2[label="fushan_040914"] du[label="data.to.use"] du2[label="data.to.use_minN1"] spdat->UU spdat->US spdat->C20S spdat->C20U spdat->mor1 [label="slow.grow.mort.onesp"] start->dat123a dat123a->dat123b dat123b->dat123 dat123->gdata gdata->mortgrow [label="(Condit,2004)"] mortgrow->mortgrow123[label="add transgrow"] mortgrow123->fu fu2->du du->du2 } ``` | | fushan.tree1-3 || fs.stem1-3 || |------|------|------|------|------| | steps|sample size|species|sample size|species| |oringinal|147450|112|264892|111| |dat123a|147450|112|170180|111| |dat123b|98758|107|96072|107| |dat123|89752|92|87550|92| |gdata|88845|92|85205|92| |mortgrow|86216|92|84113|92| |mortgrow_123|82841|91|77835|92|