nsRFA
Non-Supervised Regional Frequency Analysis

Package index
Search the nsRFA package
Vignettes
Functions
359
Source code
36
Man pages
31
Home
/
CRAN
/
nsRFA
/
demo/feh_floods02.R

demo/feh_floods02.R
In nsRFA: Non-Supervised Regional Frequency Analysis 

# ------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------- #
#                     ANALYSIS OF FEH DATA: CLASSIFICATION-VARIABLES                    #
# ------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------- #


# Data loading:

data(FEH1000)



# Criteria used in the FEH to choose stations for pooling groups:
#  n>7
#  area, saar and bfihost are known
#  urbext<0.025
#  area>0.5 \n

n <- tapply(am[,4],am[,1],length)
urbext <- cd[,"urbext1990"]
area <- cd[,"dtm_area"]
cd696 <- cd[(!is.nan(cd[,"dtm_area"]))&(!is.nan(cd[,"saar"]))&(!is.nan(cd[,"bfihost"]))&(n>7)&(urbext<0.025)&(area>0.5),]

fac <- factor(am[,"number"],levels=cd696[,"number"])
am696 <- am[!is.na(fac),]
nlevels(as.factor(am696[,"number"]))


# FEH classification variables:
#  lnAREA, lnSAAR, BFIHOST.

lcv <- tapply(am696[,4],am696[,1],LCV)
lca <- tapply(am696[,4],am696[,1],LCA)

lnAREA <- log(cd696[,"dtm_area"])
lnSAAR <- log(cd696[,"saar"])
BFIHOST <- cd696[,"bfihost"]

lcvVSclassvarFEH <- lm(lcv ~ lnAREA + lnSAAR + BFIHOST)
lcaVSclassvarFEH <- lm(lca ~ lnAREA + lnSAAR + BFIHOST)

graphics.lm <- function (regr) {
 par(mfrow=c(1,2), cex=1)
  plot(regr$fitted.values,regr$model[,1],pch=".",cex=3,xlab="Originals",ylab="Fitted")
  abline(0,1,lty=3)
  intervals <- predinterval.lm(regr)
  intervals <- intervals[order(intervals[,1]),]
  lines(intervals[,c(1,2)],lty=2)
  lines(intervals[,c(1,3)],lty=2)
  Rsq <- signif(R2.lm(regr),3)
  rmse <- signif(RMSE.lm(regr),3)
  rmsep <- signif(RMSEP(regr$model[,1],regr$fitted.values)*100,3)
  mtext(paste("R2 = ",Rsq[1]),3,-1.5,adj=0.05,cex=1)
  mtext(paste("R2adj = ",Rsq[2]),3,-2.5,adj=0.05,cex=1)
  mtext(paste("RMSE = ",rmse),1,-2.5,adj=0.95,cex=1)
  mtext(paste("RMSEP = ",rmsep,"%"),1,-1.5,adj=0.95,cex=1)

  normplot(regr$residuals,pch=".",cex=3,xlab="Residuals")
 par(mfrow=c(1,1),cex=1)
 title(main=paste(names(regr$model[1]),"~",paste(names(regr$model[-1]), collapse=", ")),cex.main=1,font.main=1)
}

graphics.lm(lcvVSclassvarFEH)
par(ask = interactive())

graphics.lm(lcaVSclassvarFEH)



# Alternative classification variables:

NGRX <- cd696[,"ihdtm_ngr_x"]
NGRY <- cd696[,"ihdtm_ngr_y"]
AREA <- cd696[,"dtm_area"]
SAAR <- cd696[,"saar"]
SPRHOST <- cd696[,"sprhost"]
FARL <- cd696[,"farl"]
RMED1D <- cd696[,"rmed_1d"]
RMED2D <- cd696[,"rmed_2d"]
RMED1H <- cd696[,"rmed_1h"]
SMDBAR <- cd696[,"smdbar"]
PROPWET <- cd696[,"propwet"]
LDP <- cd696[,"ldp"]
DPLBAR <- cd696[,"dplbar"]
ALTBAR <- cd696[,"altbar"]
DPSBAR <- cd696[,"dpsbar"]
ASPBAR <- cd696[,"aspbar"]
ASPVAR <- cd696[,"aspvar"]
URBEXT <- cd696[,"urbext1990"]

potCVnat <- data.frame(cbind(NGRX,NGRY,AREA,SAAR,BFIHOST,SPRHOST,FARL,RMED1D,RMED2D,RMED1H,
                             SMDBAR,PROPWET,LDP,DPLBAR,ALTBAR,DPSBAR,ASPBAR,ASPVAR,URBEXT))
potCVln <- log(potCVnat[-c(17,19)]); names(potCVln) <- paste("ln",names(potCVln),sep="")
potCV <- cbind(potCVnat,potCVln)

# bestregressions <- function(dip,ind) {
#  Y <- as.numeric(dip)
#  X <- ind
#  Sy <- var(Y)
#  Sx <- var(X)
#  Sxy <- var(X,Y)
#  Dm.mat <- Sx
#  Dm.H <- Sxy %*% t(Sxy)/Sy
#  require(subselect)
#  Dm.leaps <- leaps(Dm.mat, kmin=1, kmax=6, H=Dm.H, r=1, nsol=3)
#  Dm.leaps
#  for(i in 6:1) {for(j in 1:3) {print(colnames(X)[Dm.leaps$subsets[j,c(1:6),i]])}}
# }

# bestregressions(lcv,potCVnat[-19])
# bestregressions(lcv,potCVln)
# bestregressions(lcv,cbind(potCVnat[-19],potCVln[-c(1,4:5,7:13)]))
# bestregressions(lcv,cbind(potCVnat[-c(1,4:5,7:13,19)],potCVln))
# bestregressions(lcv,cbind(potCVnat[-c(1:2,19)],potCVln))
names(cbind(potCVnat[-c(1:2,19)],potCVln))

graphics.lm(lm(lcv ~ lnAREA + lnSAAR + lnRMED2D,data=potCV))
graphics.lm(lm(lcv ~ DPLBAR + lnSAAR + lnRMED2D,data=potCV))
graphics.lm(lm(lcv ~ DPLBAR + ALTBAR + lnSAAR + lnRMED2D,data=potCV))
graphics.lm(lm(lcv ~ RMED1D + ALTBAR + lnAREA + lnSAAR + lnRMED2D + lnALTBAR,data=potCV))

print(prt.lm(lm(lcv ~ DPLBAR + lnSAAR + lnRMED2D,data=potCV)))
print(prt.lm(lm(lcv ~ lnAREA + lnSAAR + lnRMED1D + lnRMED2D,data=potCV)))
print(prt.lm(lm(lcv ~ DPLBAR + ALTBAR + lnSAAR + lnRMED2D + lnALTBAR,data=potCV)))
print(prt.lm(lm(lcv ~ RMED1D + ALTBAR + lnAREA + lnSAAR + lnRMED2D + lnALTBAR,data=potCV)))


# Distance matrices and Mantel test approach: best models

# distancesAD <- AD.dist(am696[,4],am696[,1])                                   # 1 hour
# distancesCVnat <- data.frame(apply(potCVnat,2,dist))                          
# distancesCVln <- data.frame(apply(potCVln,2,dist))                            
# distancesCV <- data.frame(apply(cbind(potCVnat[-c(1:2,19)],potCVln),2,dist))  

# bestregressions(as.numeric(distancesAD),distancesCVnat)  
# bestregressions(as.numeric(distancesAD),distancesCVln)   
# bestregressions(as.numeric(distancesAD),distancesCV)     


# graphics.lm(lm(distancesAD ~ FARL + RMED2D + lnSAAR,data=distancesCV))
# graphics.lm(lm(distancesAD ~ FARL + SMDBAR + DPLBAR + lnNGRY + lnSAAR + lnRMED2D,data=distancesCV))
# graphics.lm(lm(distancesAD ~ DPLBAR + lnSAAR + lnRMED2D,data=distancesCV))


# Distance matrices and Mantel test approach: mantel test

# regr1CV <- lm(distancesAD ~ SMDBAR,data=distancesCV)
# regr2CV <- lm(distancesAD ~ RMED2D + lnSAAR,data=distancesCV)
# regr3CV <- lm(distancesAD ~ FARL + RMED2D + lnSAAR,data=distancesCV)
# regr4CV <- lm(distancesAD ~ FARL + RMED2D + DPLBAR + lnSAAR,data=distancesCV)
# regr5CV <- lm(distancesAD ~ FARL + RMED2D + lnNGRY + lnSAAR + lnDPLBAR,data=distancesCV)
# regr6CV <- lm(distancesAD ~ FARL + SMDBAR + DPLBAR + lnNGRY + lnSAAR + lnRMED2D,data=distancesCV)

# mantel.lm(regr1CV, Nperm=100)   # 2 min 
# mantel.lm(regr2CV, Nperm=100)   # 4 min
# mantel.lm(regr3CV, Nperm=100)   # 6 min
# mantel.lm(regr4CV, Nperm=100)   # 6 min
# mantel.lm(regr5CV, Nperm=100)   # 6 min
# mantel.lm(regr6CV, Nperm=100)   # 6 min
# mantel.lm(lm(distancesAD ~ RMED2D + RMED1D,data=distancesCV), Nperm=100)
# prt.lm(lm(lcv ~ FARL + SMDBAR + DPLBAR + lnNGRY + lnSAAR + lnRMED2D,data=potCV))



# Distance matrices and Mantel test approach: regression statistics

# R2.lm(regr1CV)
# R2.lm(regr2CV)
# R2.lm(regr3CV)
# R2.lm(regr4CV)
# R2.lm(regr5CV)
# R2.lm(regr6CV)

# RMSEP(regr1CV$model[,1],regr1CV$fitted.values)*100
# RMSEP(regr2CV$model[,1],regr2CV$fitted.values)*100
# RMSEP(regr3CV$model[,1],regr3CV$fitted.values)*100
# RMSEP(regr4CV$model[,1],regr4CV$fitted.values)*100
# RMSEP(regr5CV$model[,1],regr5CV$fitted.values)*100
# RMSEP(regr6CV$model[,1],regr6CV$fitted.values)*100



# Formation of groups and homogeneity.
# Classification variables:
#  0) lnAREA, lnSAAR, BFIHOST (same weights + lnAREA/sqrt(2))\n

# term1 <- log(cd696[,"dtm_area"])/(sd(log(cd696[,"dtm_area"]))*sqrt(2))
# term2 <- log(cd696[,"saar"])/sd(log(cd696[,"saar"]))
# term3 <- cd696[,"bfihost"]/sd(cd696[,"bfihost"])

# roi.cd <- data.frame(cbind(term1,term2,term3))
# row.names(roi.cd) <- cd696[,"number"]

# roi00.50year <- new.env()
# for(i in 1:696) {
# print(paste(i,"/ 696"))
# assign(as.character(row.names(roi.cd)[i]), roi.st.year(roi.cd[i,],as.data.frame(roi.cd),
#         row.names(roi.cd),am696[,"am"],am696[,"number"],test="HW and AD",station.year=250,Nsim=100), env=roi00.50year)
# }
# roi00.50year <- as.list(roi00.50year)

# estrai.region <- function (x) {x$region}
# estrai.test <- function (x) {x$test}
# regioni.50year.0 <- sapply(roi00.50year, estrai.region)
# test.50year.0 <- sapply(roi00.50year, estrai.test)
# mL.50year.0 <- mean(sapply(regioni.50year.0,length))
# mH1.50year.0 <- mean(test.50year.0["H1",])
# mH2.50year.0 <- mean(test.50year.0["H2",])
# mAD.50year.0 <- median(test.50year.0["P",])
# gH1gr2.50year.0 <- sum(test.50year.0["H1",]>2)/696
# gH1gr4.50year.0 <- sum(test.50year.0["H1",]>4)/696
# gH2gr2.50year.0 <- sum(test.50year.0["H2",]>2)/696
# gH2gr4.50year.0 <- sum(test.50year.0["H2",]>4)/696
# gADgr99.50year.0 <- sum(test.50year.0["P",]>.99)/696
# gADgr95.50year.0 <- sum(test.50year.0["P",]>.95)/696

# table.0 <- signif(c(mL.50year.0,mH1.50year.0,mH2.50year.0,mAD.50year.0,
#            gH1gr2.50year.0*100,gH1gr4.50year.0*100,gH2gr2.50year.0*100,gH2gr4.50year.0*100,
#            gADgr95.50year.0*100,gADgr99.50year.0*100),3)
# names(table.0) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
#                     "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
# print(table.0)

table.0 <- c(11.20,2.87,1.53,0.99,62.50,25.10,33.80,5.89,68.20,46.00)
names(table.0) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                   "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
print(table.0)



# Formation of groups and homogeneity. 
# Classification variables:
#  1) DPLBAR, lnSAAR, lnRMED2D (with the same weight)

# term1 <- cd696[,"dplbar"]/sd(cd696[,"dplbar"])
# term2 <- log(cd696[,"saar"])/sd(log(cd696[,"saar"]))
# term3 <- log(cd696[,"rmed_2d"])/sd(log(cd696[,"rmed_2d"]))

# roi.cd <- data.frame(cbind(term1,term2,term3))
# row.names(roi.cd) <- cd696[,"number"]

# roi.st.year(roi.cd["54088",],roi.cd,row.names(roi.cd),am696[,"am"],
#             am696[,"number"],test="HW and AD",station.year=250,Nsim=500)

# roi01.50year <- new.env()
# for(i in 1:696) {
# print(paste(i,"/ 696"))
# assign(as.character(row.names(roi.cd)[i]), roi.st.year(roi.cd[i,],as.data.frame(roi.cd),
#         row.names(roi.cd),am696[,"am"],am696[,"number"],test="HW and AD",station.year=250,Nsim=100), env=roi01.50year)
# }
# roi01.50year <- as.list(roi01.50year)

# estrai.region <- function (x) {x$region}
# estrai.test <- function (x) {x$test}
# regioni.50year.1 <- sapply(roi01.50year, estrai.region)
# test.50year.1 <- sapply(roi01.50year, estrai.test)
# mL.50year.1 <- mean(sapply(regioni.50year.1,length))
# mH1.50year.1 <- mean(test.50year.1["H1",])
# mH2.50year.1 <- mean(test.50year.1["H2",])
# mAD.50year.1 <- median(test.50year.1["P",])
# gH1gr2.50year.1 <- sum(test.50year.1["H1",]>2)/696
# gH1gr4.50year.1 <- sum(test.50year.1["H1",]>4)/696
# gH2gr2.50year.1 <- sum(test.50year.1["H2",]>2)/696
# gH2gr4.50year.1 <- sum(test.50year.1["H2",]>4)/696
# gADgr99.50year.1 <- sum(test.50year.1["P",]>.99)/696
# gADgr95.50year.1 <- sum(test.50year.1["P",]>.95)/696

# table.1 <- signif(c(mL.50year.1,mH1.50year.1,mH2.50year.1,mAD.50year.1,
#            gH1gr2.50year.1*100,gH1gr4.50year.1*100,gH2gr2.50year.1*100,gH2gr4.50year.1*100,
#            gADgr95.50year.1*100,gADgr99.50year.1*100),3)
# names(table.1) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
#                    "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
# print(table.1)
table.1 <- c(11.30,2.94,1.67,0.99,60.90,27.00,38.20,8.62,65.80,44.70)
names(table.1) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                   "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
print(table.1)





# Formation of groups and homogeneity.
# Classification variables:
#  2) DPLBAR, lnSAAR, lnRMED2D (weighted with regression coefficients)\n

coeff <- lm(lcv ~ DPLBAR + lnSAAR + lnRMED2D,data=potCV)$coefficients[-1]
print(coeff)

# term1 <- cd696[,"dplbar"]*coeff[1]
# term2 <- log(cd696[,"saar"])*coeff[2]
# term3 <- log(cd696[,"rmed_2d"])*coeff[3]

# roi.cd <- data.frame(cbind(term1,term2,term3))
# row.names(roi.cd) <- cd696[,"number"]

# roi02.50year <- new.env()
# for(i in 1:696) {
# print(paste(i,"/ 696"))
# assign(as.character(row.names(roi.cd)[i]), roi.st.year(roi.cd[i,],as.data.frame(roi.cd),
#         row.names(roi.cd),am696[,"am"],am696[,"number"],test="HW and AD",station.year=250,Nsim=100), env=roi02.50year)
# }
# roi02.50year <- as.list(roi02.50year)

# estrai.region <- function (x) {x$region}
# estrai.test <- function (x) {x$test}
# regioni.50year.2 <- sapply(roi02.50year, estrai.region)
# test.50year.2 <- sapply(roi02.50year, estrai.test)
# mL.50year.2 <- mean(sapply(regioni.50year.2,length))
# mH1.50year.2 <- mean(test.50year.2["H1",])
# mH2.50year.2 <- mean(test.50year.2["H2",])
# mAD.50year.2 <- median(test.50year.2["P",])
# gH1gr2.50year.2 <- sum(test.50year.2["H1",]>2)/696
# gH1gr4.50year.2 <- sum(test.50year.2["H1",]>4)/696
# gH2gr2.50year.2 <- sum(test.50year.2["H2",]>2)/696
# gH2gr4.50year.2 <- sum(test.50year.2["H2",]>4)/696
# gADgr99.50year.2 <- sum(test.50year.2["P",]>.99)/696
# gADgr95.50year.2 <- sum(test.50year.2["P",]>.95)/696

# table.2 <- signif(c(mL.50year.2,mH1.50year.2,mH2.50year.2,mAD.50year.2,
#            gH1gr2.50year.2*100,gH1gr4.50year.2*100,gH2gr2.50year.2*100,gH2gr4.50year.2*100,
#            gADgr95.50year.2*100,gADgr99.50year.2*100),3)
# names(table.2) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
#                     "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
# print(table.2)

table.2 <- c(11.20,2.95,1.64,0.99,57.20,28.60,37.90,9.34,66.40,45.10) # weigths doesnt change the result
names(table.2) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                   "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
print(table.2)


# Formation of groups and homogeneity.
# Classification variables:
#  3) FARL, RMED2D, lnSAAR (same weight)

# term1 <- cd696[,"farl"]/sd(cd696[,"farl"])
# term2 <- cd696[,"rmed_2d"]/sd(cd696[,"rmed_2d"])
# term3 <- log(cd696[,"saar"])/sd(log(cd696[,"saar"]))

# roi.cd <- data.frame(cbind(term1,term2,term3))
# row.names(roi.cd) <- cd696[,"number"]

# roi03.50year <- new.env()
# for(i in 1:696) {
# print(paste(i,"/ 696"))
# assign(as.character(row.names(roi.cd)[i]), roi.st.year(roi.cd[i,],as.data.frame(roi.cd),
#         row.names(roi.cd),am696[,"am"],am696[,"number"],test="HW and AD",station.year=250,Nsim=100), env=roi03.50year)
# }
# roi03.50year <- as.list(roi03.50year)

# estrai.region <- function (x) {x$region}
# estrai.test <- function (x) {x$test}
# regioni.50year.3 <- sapply(roi03.50year, estrai.region)
# test.50year.3 <- sapply(roi03.50year, estrai.test)
# mL.50year.3 <- mean(sapply(regioni.50year.3,length))
# mH1.50year.3 <- mean(test.50year.3["H1",])
# mH2.50year.3 <- mean(test.50year.3["H2",])
# mAD.50year.3 <- median(test.50year.3["P",])
# gH1gr2.50year.3 <- sum(test.50year.3["H1",]>2)/696
# gH1gr4.50year.3 <- sum(test.50year.3["H1",]>4)/696
# gH2gr2.50year.3 <- sum(test.50year.3["H2",]>2)/696
# gH2gr4.50year.3 <- sum(test.50year.3["H2",]>4)/696
# gADgr99.50year.3 <- sum(test.50year.3["P",]>.99)/696
# gADgr95.50year.3 <- sum(test.50year.3["P",]>.95)/696

# table.3 <- signif(c(mL.50year.3,mH1.50year.3,mH2.50year.3,mAD.50year.3,
#            gH1gr2.50year.3*100,gH1gr4.50year.3*100,gH2gr2.50year.3*100,gH2gr4.50year.3*100,
#            gADgr95.50year.3*100,gADgr99.50year.3*100),3)
# names(table.3) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
#                     "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
# print(table.3)

table.3 <- c(11.00,3.37,1.60,0.99,67.20,35.10,34.60,8.05,67.40,46.30)
names(table.3) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                   "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
print(table.3)




# Formation of groups and homogeneity.
# Classification variables:
#  0) lnAREA, lnSAAR, BFIHOST (same weight + lnAREA/sqrt(2))
#  1) DPLBAR, lnSAAR, lnRMED2D (with the same weight)
#  2) DPLBAR, lnSAAR, lnRMED2D (weighted with regression coefficients)
#  3) FARL, RMED2D, lnSAAR (with the same weight)\n

table.0123 <- rbind(table.0,table.1,table.2,table.3)
print(table.0123)


# Formation of groups and homogeneity (other cases).
# Classification variables:

# Classification variables:
#  0) lnAREA, lnSAAR, BFIHOST (same weight + lnAREA/sqrt(2))
#  1) DPLBAR, lnSAAR, lnRMED2D (with the same weight)
#  2) DPLBAR, lnSAAR, lnRMED2D (weighted with regression coefficients)
#  3) FARL, RMED2D, lnSAAR (with the same weight)
#  4) FARL, SMDBAR, DPLBAR, lnNGRY, lnSAAR, lnRMED2D (with the same weight)
#  5) RMED1D, ALTBAR, lnAREA, lnSAAR, lnRMED2D, lnALTBAR (with the same weight)\n

poolingroupsum <- function(roi.cd) {
 roi.50year <- new.env()
 for(i in 1:696) {
 print(paste(i,"/ 696"))
 assign(as.character(row.names(roi.cd)[i]), roi.st.year(roi.cd[i,],as.data.frame(roi.cd),
         row.names(roi.cd),am696[,"am"],am696[,"number"],test="HW and AD",station.year=250,Nsim=100), env=roi.50year)
 }
 roi.50year <- as.list(roi.50year)

 estrai.region <- function (x) {x$region}
 estrai.test <- function (x) {x$test}
 regioni.50year <- sapply(roi.50year, estrai.region)
 test.50year <- sapply(roi.50year, estrai.test)
 mL.50year <- mean(sapply(regioni.50year,length))
 mH1.50year <- mean(test.50year["H1",])
 mH2.50year <- mean(test.50year["H2",])
 mAD.50year <- median(test.50year["P",])
 gH1gr2.50year <- sum(test.50year["H1",]>2)/696
 gH1gr4.50year <- sum(test.50year["H1",]>4)/696
 gH2gr2.50year <- sum(test.50year["H2",]>2)/696
 gH2gr4.50year <- sum(test.50year["H2",]>4)/696
 gADgr99.50year <- sum(test.50year["P",]>.99)/696
 gADgr95.50year <- sum(test.50year["P",]>.95)/696
 
 table <- signif(c(mL.50year,mH1.50year,mH2.50year,mAD.50year,
            gH1gr2.50year*100,gH1gr4.50year*100,gH2gr2.50year*100,gH2gr4.50year*100,
            gADgr95.50year*100,gADgr99.50year*100),3)
 names(table) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                     "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
 table
}

termFARL <- cd696[,"farl"]/sd(cd696[,"farl"])
termSMDBAR <- cd696[,"smdbar"]/sd(cd696[,"smdbar"])
termDPLBAR <- cd696[,"dplbar"]/sd(cd696[,"dplbar"])
termlnNGRY <- log(cd696[,"ihdtm_ngr_y"])/sd(log(cd696[,"ihdtm_ngr_y"]))
termlnSAAR <- log(cd696[,"saar"])/sd(log(cd696[,"saar"]))
termlnRMED2D <- log(cd696[,"rmed_2d"])/sd(log(cd696[,"rmed_2d"]))
termRMED1D <- cd696[,"rmed_1d"]/sd(cd696[,"rmed_1d"])
termALTBAR <- cd696[,"altbar"]/sd(cd696[,"altbar"])
termlnAREA <- log(cd696[,"dtm_area"])/sd(log(cd696[,"dtm_area"]))
termlnALTBAR <- log(cd696[,"altbar"])/sd(log(cd696[,"altbar"]))

# table.4 <- poolingroupsum(data.frame(cbind(termFARL,termSMDBAR,termDPLBAR,termlnNGRY,termlnSAAR,termlnRMED2D),
#            row.names=cd696[,"number"]))
# table.5 <- poolingroupsum(data.frame(cbind(termRMED1D,termALTBAR,termlnAREA,termlnSAAR,termlnRMED2D,termlnALTBAR),
#            row.names=cd696[,"number"]))

table.4 <- c(11.10,2.68,1.56,0.98,61.20,24.10,35.20,6.03,61.90,40.40)
names(table.4) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                   "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")
table.5 <- c(11.1,2.91,1.77,0.99,62.5,24.9,39.9,9.05,65.9,44.8)
names(table.4) <- c("Avg. n sites","m(H1)","m(H2)","med(p(AD))","% H1>2","% H1>4",
                   "% H2>2","% H2>4","% p(AD)>0.95","% p(AD)>0.99")

print(rbind(table.4,table.5))
print(rbind(table.0123,table.4,table.5))

# ------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------- #
#                                       THE END                                         #
# ------------------------------------------------------------------------------------- #
# ------------------------------------------------------------------------------------- #

Try the nsRFA package in your browser

Any scripts or data that you put into this service are public.

nsRFA documentation built on May 29, 2024, 1:14 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close