################################ # # Various analysis for P. superfusa data # # - Questions about patterns of survivorship # - Descriptions of growth / shrinkage # # - Various items of data handling and summary statistics # # *** Code could have commenting improved (12 June 2015) # ################################ ### # Statistical To-Do list # # 5 June 2015 # # &x&x&x Calculate number of colonies that disappeared (size = 0) and then reappeared xxx DONE # - This is a precursor to the chapter on cryptic tissues # - Operational limits of photographic methods # # &x&x&x Partition the growth from 'recruitment' into colony in a new location vs colony in a place with a recent adult death (i.e., a #3, as described in 'with_zero') # # &x&x&x ** Check how many 'mortality' events in logistic regression are included because of these colonies ** # - We don't care about guys that may re-grow later from 'death'. We are defining them operatively as mortality events. # # &x&x&x Re-make graphs by site # # # &x&x&x Test for model assumptions of ANCOVA -- residual plots # # &x&x&x Repeat ANCOVA using colonies only once # # &x&x&x Test for effects of site as covariate in ANCOVA xxx DONE -- Site is an independent effect (model5) # # &x&x&x Perform posthoc to look for which years and which sites are different # # &x&x&x Double-check that we set all values <1 to be equal to zero ### ###### # CORE ANALYSIS -- size-dependence of survivorship ###### ################################ # Use of logistic regression for P. superfusa data ################################ superfusa = read.csv("PsuperfusaRAW_forR.csv", header = TRUE) # Flexibility here to decide on the time horizon for survivorship # - Select 'start' and 'end' time for consideration of survivorship time_start = 5 # time point of the first reference -- 2009="5", 2010="6", 2011="7", 2012="8" time_start = 8 # time point of the last reference -- 2009="5", 2010="6", 2011="7", 2012="8" n = dim(superfusa)[1] # Coding for survivorship options across defined time window # live = 1 / die = 0 Colony was above min_size at 'start' and at 'end' ("survivor") # live = 0 / die = 1 Colony was above min_size at 'start' and below min_size at 'end' ("dead") # live = 0 / die = 0 Colony was below min_size at 'start' and at 'end' ("not present or considered in this time window") survivors = array(NA,c(n,3),dimnames=list(c(1:n),c("size","live","die"))) # "size" defined as size at time_start # Minimum size indicating that any colony reported as smaller than this is considered not reliable data and reported as 0-size min_size = 1 for (i in 1:n) { survivors[i,1] = superfusa[i,time_start] if(is.na(sum(superfusa[i,5:8]))==FALSE) { if(superfusa[i,time_start] >= min_size) { if(superfusa[i,time_end] >= min_size) { survivors[i,2] = 1 survivors[i,3] = 0 } else { survivors[i,2] = 0 survivors[i,3] = 1 } } else { # Analysis focuses on mortality, and thus new recruits are not involved statistically survivors[i,2] = 0 survivors[i,3] = 0v } } else { survivors[i,2] = 0 survivors[i,3] = 0 } } # Creating a one-dimensional vector with colony-specific codes of survival # 1 = survival over time window # 0 = death over time window # NA = not considered over time window (i.e., new recruit or not present in either time period) temp_surv = array(NA, n) for (i in 1:n) { if (survivors[i,2] + survivors[i,3] > 0) { temp_surv[i] = survivors[i,2] } } # Visualizing survival fate as a function of starting size plot(temp_surv ~ survivors[,1]) # Logistic regression model (using 'binomial' link function) # *** Two versions of data entry tried in a effort to better understand function glm() # *** - Should produce identical statistical output # *** - First version can handle non-binary input data (e.g., each replicate could have multiple binary trials reported) # *** - Second version used only binary input data (i.e., each datum is '0' or '1' as a single coin flip) y = cbind(survivors[,2], survivors[,3]) model = glm(y ~ survivors[,1], binomial) y2 = temp_surv model2 = glm(y2 ~ survivors[,1], binomial) summary(model) summary(model2) # Graphical output -- best-fit model parameters and raw data across size model_intercept = model$coefficients[1] model_slope = model$coefficients[2] max_size = max(superfusa[,time_start], na.rm=TRUE) size = seq(0,max_size,by=1) log_odds_ratio = model_slope*size + model_intercept prob_s = exp(log_odds_ratio) / (1 + exp(log_odds_ratio)) plot(size, prob_s, ylim=c(0,1)) points(survivors[,2] ~ survivors[,1], col="red") #################### # END logistic regression of size vs survivorship #################### ### # Relationship between changes across years ### n = dim(superfusa)[1] # 'deltaXXYY' arrays report the proportional change in area of survivors and non-recruits # - reported as the difference of ln( area_tZ ) values, which is equivalent to ln( area_t1 / area_t0 ) # - presents patterns of proportional area change (e.g., colony 1 reduced to 45% of original size [or more accurately ln(0.45) of original size]) delta0910 = array(NA,n) delta1011 = array(NA,n) delta1112 = array(NA,n) # 'abs_deltaXXYY' arrays report the absolute change in area of all colonies, independent of fate # - reported as the difference of area in units of cm^2 # * recruits are reported as the size of the recruit (which is the size of the recruit minus 0, since nothing was there before) # * dying colonies are reported as the negative of the size of the original colony (0 minus size of original colony) # Report changes in area (absolute delta) as well as the 'fate' associated (t is time and x is size) # fate = 1 True recruit [t-1, x-1=0; t0, x0=0; t1, x1>0] # fate = 2 Resurrected recruit [t-1, x-1>0; t0, x0=0; t1, x1>0] # fate = 3 Growth [t-1, x-1>=0; t0, x0>0; t1, x1>x0] # fate = 4 Shrink [t-1, x-1>=0; t0, x0>0; t1, x1=0; t0, x0>0; t1, x1=0] abs_delta0910 = array(NA,c(n,2),dimnames=list(c(1:n),c("change_in_size(cm^2)","fate"))) abs_delta1011 = array(NA,c(n,2),dimnames=list(c(1:n),c("change_in_size(cm^2)","fate"))) abs_delta1112 = array(NA,c(n,2),dimnames=list(c(1:n),c("change_in_size(cm^2)","fate"))) # Minimum size indicating that any colony reported as smaller than this is considered not reliable data and reported as 0-size min_size = 1 for (i in 1:n) { if(is.na(sum(superfusa[i,5:8]))==FALSE) { # For time period 2009-2010 if(superfusa[i,5] > min_size && superfusa[i,6] > min_size) { delta0910[i] = log(superfusa[i,6])-log(superfusa[i,5]) abs_temp = superfusa[i,6]-superfusa[i,5] # Fate 3 if(abs_temp >= 0) abs_delta0910[i,] = c(abs_temp,3) # Fate 4 if(abs_temp < 0) abs_delta0910[i,] = c(abs_temp,4) } if(superfusa[i,5] > min_size && superfusa[i,6] < min_size) { abs_temp = (-1)*superfusa[i,5] # Fate 5 abs_delta0910[i,] = c(abs_temp,5) } if(superfusa[i,5] < min_size && superfusa[i,6] > min_size) { abs_temp = superfusa[i,6] # Fate 1 (no information in first pair of years to discriminate between 'true' and 'resurrected' recruits) abs_delta0910[i,] = c(abs_temp,1) } # For time period 2010-2011 if(superfusa[i,6] > min_size && superfusa[i,7] > min_size) { delta1011[i] = log(superfusa[i,7])-log(superfusa[i,6]) abs_temp = superfusa[i,7]-superfusa[i,6] # Fate 3 if(abs_temp >= 0) abs_delta1011[i,] = c(abs_temp,3) # Fate 4 if(abs_temp < 0) abs_delta1011[i,] = c(abs_temp,4) } if(superfusa[i,6] > min_size && superfusa[i,7] < min_size) { abs_temp = (-1)*superfusa[i,6] # Fate 5 abs_delta1011[i,] = c(abs_temp,5) } if(superfusa[i,6] < min_size && superfusa[i,7] > min_size) { abs_temp = superfusa[i,7] # Fate 1 if(superfusa[i,5] < min_size) abs_delta1011[i,] = c(abs_temp,1) # Fate 2 if(superfusa[i,5] >= min_size) abs_delta1011[i,] = c(abs_temp,2) } # For time period 2011-2012 if(superfusa[i,7] > min_size && superfusa[i,8] > min_size) { delta1112[i] = log(superfusa[i,8])-log(superfusa[i,7]) abs_temp = superfusa[i,8]-superfusa[i,7] # Fate 3 if(abs_temp >= 0) abs_delta1112[i,] = c(abs_temp,3) # Fate 4 if(abs_temp < 0) abs_delta1112[i,] = c(abs_temp,4) } if(superfusa[i,7] > min_size && superfusa[i,8] < min_size) { abs_temp = (-1)*superfusa[i,7] # Fate 5 abs_delta1112[i,] = c(abs_temp,5) } if(superfusa[i,7] < min_size && superfusa[i,8] > min_size) { abs_temp = superfusa[i,8] # Fate 1 if(superfusa[i,6] < min_size) abs_delta1112[i,] = c(abs_temp,1) # Fate 2 if(superfusa[i,6] >= min_size) abs_delta1112[i,] = c(abs_temp,2) } } } # Summarize the change in absolute area (already separated by time and by fate) by site # - Provides a summed areal change by site sites = levels(superfusa[,2]) n_sites = length(sites) abs_delta_by_fate_0910 = array(NA, c(n_sites,6), dimnames=list(c(sites),c("True_recruit","Resurrected_recruit","Growth","Shrink","Death","Net_change"))) abs_delta_by_fate_1011 = array(NA, c(n_sites,6), dimnames=list(c(sites),c("True_recruit","Resurrected_recruit","Growth","Shrink","Death","Net_change"))) abs_delta_by_fate_1112 = array(NA, c(n_sites,6), dimnames=list(c(sites),c("True_recruit","Resurrected_recruit","Growth","Shrink","Death","Net_change"))) for(i in 1:n_sites) { abs_delta_by_fate_0910[i,1] = tapply(abs_delta0910[which(superfusa[,2]==sites[i]),1],abs_delta0910[which(superfusa[,2]==sites[i]),2],sum)[1] abs_delta_by_fate_0910[i,c(3:5)] = tapply(abs_delta0910[which(superfusa[,2]==sites[i]),1],abs_delta0910[which(superfusa[,2]==sites[i]),2],sum)[c(2:4)] abs_delta_by_fate_1011[i,c(1:5)] = tapply(abs_delta1011[which(superfusa[,2]==sites[i]),1],abs_delta1011[which(superfusa[,2]==sites[i]),2],sum) abs_delta_by_fate_1112[i,c(1:5)] = tapply(abs_delta1112[which(superfusa[,2]==sites[i]),1],abs_delta1112[which(superfusa[,2]==sites[i]),2],sum) abs_delta_by_fate_0910[i,6] = sum(abs_delta_by_fate_0910[i,c(1,3:5)]) abs_delta_by_fate_1011[i,6] = sum(abs_delta_by_fate_1011[i,c(1:5)]) abs_delta_by_fate_1112[i,6] = sum(abs_delta_by_fate_1112[i,c(1:5)]) } # Summaries at a per-quad level (dividing by number of quads per site [n=10]) # - Recall that each quadrat covers an area of 60cm x 90cm, or 5400cm^2 # * Thus, proportional areal estimates of change need to account for the per-quad areal coverage #D = 1 # Reports the total areal change per site in all quads D = 10 # Reports the average areal change per quadrat for each site #D = 10 * 5400 # Reports the proportional areal change per quadrat for each site abs_delta_by_fate_0910/D abs_delta_by_fate_1011/D abs_delta_by_fate_1112/D #################### # END relationship of changes among years #################### ### # Statistical test of pattern of growth # - Testing for autocorrelation in growth rate; if you grew quickly in the last time window, is that linked to current time window? # * testing for site effects (though results show little effect of 'site' on relationship) # - There are two tests here, given that there are data from four years (Year1-->Year2 vs Year2-->Year3 & Year2-->Year3 vs Year3-->Year4) ### # ANCOVA model 1 -- testing for an interaction (i.e., different slopes for different sites) y=aov(delta1011 ~ delta0910 * superfusa$Site) summary(y) # ANCOVA model 2 -- If no significant interaction, test for group-specific differences in intercept y=aov(delta1011 ~ delta0910 + superfusa$Site) summary(y) # ANCOVA model 3 -- If no significant effect of Site after accounting for covariate, proceed to a typical linear regression y=aov(delta1011 ~ delta0910) summary(y) # or...use the lm() function, which gives you the regression results that you want <> y=lm(delta1011~delta0910) summary(y) # Graphical view of covariate effect (without site grouping) plot(delta1011~delta0910) abline(y) points(0,0, col="red", pch=19) # Test of residuals, etc plot(y) ### # Relationship of 2012/2011 transition vs 2011/2010 transition plot(delta1112~delta1011) # ANCOVA model 1 -- testing for an interaction (i.e., different slopes for different sites) y=aov(delta1112 ~ delta1011 * superfusa$Site) summary(y) # ANCOVA model 2 -- If no significant interaction, test for group-specific differences in intercept y=aov(delta1112 ~ delta1011 + superfusa$Site) summary(y) # ANCOVA model 3 -- If no significant effect of Site after accounting for covariate, proceed to a typical linear regression y=aov(delta1112 ~ delta1011) summary(y) # or...use the lm() function, which gives you the regression results that you want <> y=lm(delta1112~delta1011) summary(y) # Graphical view of covariate effect (without site grouping) plot(delta1112~delta1011) abline(y) points(0,0, col="red", pch=19) # Test of residuals, etc plot(y) ## Plot the year-change comparisons on the same graph plot(delta1011~delta0910, pch=16, col="pink", xlim=c(-4,3), ylim=c(-3,3), xlab="Proportional growth, window 1", ylab="Proportional growth, window 2") points(delta1112~delta1011, pch=17, col="turquoise") points(0,0, col="red", pch=19) y=lm(delta1011~delta0910) abline(y, col="pink") y=lm(delta1112~delta1011) abline(y, col="turquoise") legend(-4, 3.2, c("Growth0910 vs Growth1011", "Growth1011 vs Growth1112"), col = c("pink", "turquoise"), pch = c(16, 17)) ### # END relationship among sequential proportional growth rates ### ###### # CORE ANALYSIS -- size-dependence of growth ###### ### # Relationships between proportional growth and starting size -- are there size effects on growth? ### # Graphical exploration of the changes plot(delta0910 ~ log(superfusa[,5])) y0910=lm(delta0910 ~ log(superfusa[,5])) summary(y0910) plot(delta1011 ~ log(superfusa[,6])) y1011=lm(delta1011 ~ log(superfusa[,6])) summary(y1011) plot(delta1112 ~ log(superfusa[,7])) y1112=lm(delta1112 ~ log(superfusa[,7])) summary(y1112) # Plot all of the points on one graph plot(delta0910 ~ log(superfusa[,5]), xlim=c(0,6.5), ylim=c(-4,3)) abline(y0910) points(delta1011 ~ log(superfusa[,6]), col="red") abline(y1011, col="red") points(delta1112 ~ log(superfusa[,7]), col="blue") abline(y1112, col="blue") # Plotting all data, coding by site and year # * Version #1 -- Color reflects year site_levels=c("FR3","FR5","FR7","FR9") site_labels=c("3","5","7","9") plot(0 ~ 0, pch=".", xlim=c(0,6.5), ylim=c(-4,3), xlab="ln(starting size)", ylab="ln(size_ratio)") size=0.65 for (i in 1:4) { points(delta0910[which(superfusa$Site==site_levels[i])] ~ log(superfusa[which(superfusa$Site==site_levels[i]),5]), pch=site_labels[i], cex=size) points(delta1011[which(superfusa$Site==site_levels[i])] ~ log(superfusa[which(superfusa$Site==site_levels[i]),6]), pch=site_labels[i], col="red", cex=size) points(delta1112[which(superfusa$Site==site_levels[i])] ~ log(superfusa[which(superfusa$Site==site_levels[i]),7]), pch=site_labels[i], col="blue", cex=size) } # * Version #2 -- Color reflects site plot(0 ~ 0, pch=".", xlim=c(0,6.5), ylim=c(-4,3), xlab="ln(starting size)", ylab="ln(size_ratio)") site_colors=c(10,35,85,84) for (i in 1:4) { points(0,-i,col=site_colors[i]) points(delta0910[which(superfusa$Site==site_levels[i])] ~ log(superfusa[which(superfusa$Site==site_levels[i]),5]), pch="1", col=site_colors[i], cex=size) points(delta1011[which(superfusa$Site==site_levels[i])] ~ log(superfusa[which(superfusa$Site==site_levels[i]),6]), pch="2", col=site_colors[i], cex=size) points(delta1112[which(superfusa$Site==site_levels[i])] ~ log(superfusa[which(superfusa$Site==site_levels[i]),7]), pch="3", col=site_colors[i], cex=size) } # Conduct an ANCOVA with year-range as grouping factor and size at time 't0' as the continuous covariate # - First need to concatenate data into one file years = c(rep("0910",1151),rep("1011",1151),rep("1112",1151)) size_start = c(log(superfusa[,5]),log(superfusa[,6]),log(superfusa[,7])) size_delta = c(delta0910, delta1011, delta1112) site = as.factor(c(rep(superfusa$Site,3))) # Model 2 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with all interaction terms model2=aov(size_delta ~ as.factor(years) * size_start * site) summary(model2) # Model 3 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with all interaction terms model3=aov(size_delta ~ site + as.factor(years) * size_start) summary(model3) # Model 4 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with interaction between 'site' and 'transition year' model4=aov(size_delta ~ site * as.factor(years) + size_start) summary(model4) # Model 5 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with interaction between 'site' and 'starting size' model5=aov(size_delta ~ site * size_start + as.factor(years)) summary(model5) # Model 6 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with no interactions model6=aov(size_delta ~ as.factor(years) + site + size_start) summary(model6) TukeyHSD(model6, "as.factor(years)") TukeyHSD(model6, "site") ###### # BEGIN counting routine ###### # Test for data volume for deltas. How many colonies were 'repeat customers'? number_obs=array(NA,length(delta0910)) # In same loop, count which colonies have a non-zero value, then a zero value, then a non-zero value (photographic 'zombie') with_zero=array(NA,length(delta0910)) for (i in 1:length(delta0910)) { temp=0 if(is.na(delta0910[i])==FALSE) temp=temp+1 if(is.na(delta1011[i])==FALSE) temp=temp+1 if(is.na(delta1112[i])==FALSE) temp=temp+1 number_obs[i] = temp temp_zero = 0 # Keeps code for fates -- 0 = zero size, 1 = non-zero size, 2 = zero size after non-zero size, 3 = non-zero size after being '2', 4 = NA in the data series for (j in 5:8) { temp_size = superfusa[i,j] if (is.na(temp_size) == TRUE) { temp_zero = 4 break } if(temp_zero == 2 && temp_size > 0) temp_zero = 3 if(temp_zero == 1 && temp_size == 0) temp_zero = 2 if(temp_zero == 0 && temp_size > 0) temp_zero = 1 } with_zero[i] = temp_zero } table(number_obs) # reports the distribution of colonies with particular number of transitions table(with_zero) # reports the number of corals that existed [1], that died [2], and that 'died' and recovered [3] ###### # END counting routine ###### ##### # # Creating a new vector of delta data where each colony is represented a maximum of one time # - No non-independence among the delta data # * 'number_obs' tells us the number of delta observations for each unique coral ID # - This code will create just the subset delta data (named 'size_delta1') which can be run identically through ANCOVA scripts # ##### # Note that this will provide different quantitatively results each time you run it, given that it is a subsampling routine size_delta1=size_delta # Based on 'number_obs', we only have to modify colonies with a '2' or '3' value n_colonies = length(number_obs) for (i in 1:n_colonies) { if (number_obs[i] == 2 || number_obs[i] == 3) { temp_rows = c(i,(i+n_colonies),(i+n_colonies+n_colonies)) temp_deltas = size_delta1[temp_rows] temp_position = 4 while(is.na(temp_deltas[temp_position])==TRUE) { temp_position = sample(1:3,1) } final_deltas=c(NA,NA,NA) final_deltas[temp_position] = temp_deltas[temp_position] size_delta1[temp_rows] = final_deltas } } #mean(size_delta1,na.rm=TRUE) # Model 2 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with all interaction terms model2a=aov(size_delta1 ~ as.factor(years) * size_start * site) summary(model2a) # Model 3 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with all interaction terms model3a=aov(size_delta1 ~ site + as.factor(years) * size_start) summary(model3a) # Model 4 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with interaction between 'site' and 'transition year' model4a=aov(size_delta1 ~ site * as.factor(years) + size_start) summary(model4a) # Model 5 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with interaction between 'site' and 'starting size' model5a=aov(size_delta1 ~ site * size_start + as.factor(years)) summary(model5a) # Model 6 # Two-factor ANCOVA -- ln(prop. growth) as a function of transition year & site (grouping factors) # * starting size as covariate # * run with no interactions model6a=aov(size_delta1 ~ as.factor(years) + site + size_start) summary(model6a) TukeyHSD(model6a, "as.factor(years)") TukeyHSD(model6a, "site") # END subsampling routine ###### # # .___, # ___('v')___ # `"-\._./-"' # ^ ^ # # Kate flies away with code before # 12 June 2015 # ##### ################# # Testing relationship of adding a constant area to a colony independent of starting size # *** Need to include a 'death' term that allows colony growth to go negative # *** (Failed) attempt at this is below ################# size1=seq(1,100,by=1) y=2 size2=size1+y plot(log(size2/size1)~log(size1)) ####### # # Simulation of coral growth / death # # - Growth modeled as constant radial expansion (independent of size); g = radial growth per unit time # - Death modeled as constant proportional loss (independent of size); d = proportion of area lost per unit time # # *** Did not really work... *** # 14 May 2015 ####### g = 1 d = 0.1 r = 1 # Number of time points for growth; length of growth vector n = 100 # Vector of area values through time A = array(NA,n) # Vector of radii through time for(i in 2:n) r=c(r,(r[i-1]+g)) A[1] = pi * r[1]^2 for(i in 2:n) A[i]= ((1-d)*A[i-1]) + (pi*((r[i]^2)-(r[i-1]^2))) plot(A) ### # DOS pasture -- Old code that we don't use # 12 June 2015 ### ### # _,--''--,_ _,-'~~'-, # ( '' {}{} # }{ {}{} # }{| {}{} # }{}( ) ,__ , \ {}{} # {}{) / \' ) '\ / ~/\ (\/) # {}{/ / \ | | /| / \ }{ # }{( ) ( ) ( ) ( ) | |} # {}{\\ \\ (/ || | u|{ # }{}{\\ \\ || || \ | # |\ |\ |\ |\ (\o # -- '' -'-' ~` "" ~` ~`' ~` '"' "' -- '' -'-' ### # These are the counts of value where (i) there are no NAs for the colony ### and (ii) there were non-zero values for each of the two time points in question sum(!is.na(delta0910)) sum(!is.na(delta1011)) sum(!is.na(delta1112)) ### # Relationship of 2011/2010 transition vs 2010/2009 transition plot(delta1011~delta0910) # Model 1 # One-factor ANCOVA -- ln(prop. growth) as a function of transition year (grouping factor) # * starting size as covariate model1=aov(size_delta ~ as.factor(years) * size_start) #model1_res=model1$res ### Need to test the residuals for assumptions. ### Need to account for sites (nested design) ######### # Simulation to explore behavior of logistic regression ######### model_intercept = -0.517876 model_slope = 0.018404 size = seq(0,400,by=1) log_odds_ratio = model_slope*size + model_intercept plot(size, log_odds_ratio) prob_s = exp(log_odds_ratio) / (1 + exp(log_odds_ratio)) plot(size, prob_s, ylim=c(0,1)) ## Logistic regression n=100 model_intercept = 2 model_slope = -0.01 survivors = array(NA,c(n,3),dimnames=list(c(1:n),c("size","live","die"))) for (i in 1:n) { survivors[i,1] = runif(1,0,400) temp1 = model_intercept + (model_slope * survivors[i,1]) thresh_prob = exp(temp1) / (1 + exp(temp1)) temp_rand = runif(1) if (temp_rand < thresh_prob && i < 80) { survivors[i,2] = 1 survivors[i,3] = 0 } if (temp_rand > thresh_prob && i < 80) { survivors[i,2] = 0 survivors[i,3] = 1 } if (i >= 80) { survivors[i,2] = 0 survivors[i,3] = 0 } } plot(survivors[,2] ~ survivors[,1]) # Logistic regression model (using 'binomial' link function) y = cbind(survivors[,2], survivors[,3]) model = glm(y ~ survivors[,1], binomial) summary(model)