# This function estimates parameters assuming no inputs into the population # No missing years allowed in N; Interpolate if there are missing years DHMethod_function(N, filtlen=4) #default filter length is 4 { ####Check if there is enough data; break if not if(length(N)<(filtlen+8)) { print("Error: not enough data"); break } ####Set up running sum filter and runsum runsumfilt _ rep(1,filtlen) runsum _ filter(N,runsumfilt,method="con",sides=1); runsum_runsum[is.number(runsum)] #filter func puts in trailing NAs; remove runsum_as.vector(runsum) #remove Splus timeseries data format if(is.element(0,runsum)) {print("Error: Zeros in runsum"); break} ####Get logratios of N and runsum at tau=1 runlen_length(runsum) Rlogratios _ log(runsum[2:runlen])-log(runsum[1:(runlen-1)]); lenN_length(N) Nlogratios_log(N[2:lenN]/N[1:(lenN-1)]) Nlogratiosyrs_(1:length(Nlogratios))[is.finite(Nlogratios)] #indices where Nlogratios is not zero; need for trend test Nlogratios_Nlogratios[is.finite(Nlogratios)] #chuck cases where N is zero ####Do diagnostics on data #check for normality normfit_ks.gof(Nlogratios,dist="norm") print(paste("p.value for normality goodness of fit test:", format(normfit\$p.value,digits=4))) #Test for trends in data trendstat_ls.print(lsfit(Nlogratiosyrs,Nlogratios),print.it=F) trendpvalue_trendstat\$summary[7] names(trendpvalue)_NULL print(paste("p.value for trend test:",format(trendpvalue,digits=4))) #check for 1st order autocorrelation; detrend data before doing this lenNlogs_length(Nlogratios) logratiodetrend_ Nlogratios-(trendstat\$coef.table[1,1]+trendstat\$coef.table[2,1]* (1:lenNlogs)) tmpacf_cor.test(logratiodetrend[1:(lenNlogs-1)], logratiodetrend[2:lenNlogs],method="spearman") print(paste("p.value for autocorrelation test:", format(tmpacf\$p.value,digits=4))) #check for and remove outliers in runsum logratios which is used to calc muR mulsfit_lsfit(rep(1,(runlen-1)),Rlogratios,intercept=F) mulsdiag_ls.diag(mulsfit) outliercutoff_qt(.9,runlen-1) #this is a generous; qt(.975) is standard mudfits_abs(mulsdiag\$dfits) .01) print("Sigma2slp less than 0 and not zero"); if (sigma2slp < 0 & lmsigma2sum\$sigma < .01) { sigma2slp_0.00001; print("Replaced neg sigma2slp with 1E-5") } #Get sigma2np sigma2np_(var(Nlogratios)-sigma2slp)/2 #Get sigma2muR sigma2muR_(2*sigma2np/filtlen + (lenN-filtlen)*sigma2slp)/(lenN-filtlen)^2 return(list(muR=muR, sigma2slp=sigma2slp, sigma2np=sigma2np, sigma2muR=sigma2muR, n=lenN, L=filtlen)) }