Først så prøver jeg å gjenskape fjellindeksen slik den er på nettsiden. Indeksen skal være 0.56.
Fyll inn ditt eget passord og brukernavn
myUser <- "anders.kolstad@nina.no"
myPwd <- "" # hemmelig passord
Importerer data fra NI-databasen og lagrer datasettet på server
fjelldat <- NIcalc::importDatasetApi(
username = myUser,
password = myPwd,
eco = "Fjell",
indic = NULL,
year = c("1990","2000","2010","2014","2019"))
saveRDS(fjelldat, "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-1990-2019.rds")
rm(fjelldat)
rm(myPwd)
Laster tilbake datasettet:
#fjelldat <- readRDS("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-1990-2019.rds")
fjelldat <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-1990-2019.rds")
Liste over arter som inngår
(inds <- fjelldat$indicators$name)
## [1] "Heilo" "Boltit"
## [3] "Snøspurv" "Steinskvett"
## [5] "Heipiplerke" "Fjellerke"
## [7] "Ravn" "Blåstrupe"
## [9] "Ringtrost" "Lappspurv"
## [11] "Fjellrype" "Villrein"
## [13] "Lirype" "smågnagere - fjellbestander"
## [15] "Tamrein" "Vier alpint belte"
## [17] "Fjellfiltmose" "Issoleie"
## [19] "Praktdraugmose" "Fjellvalmue"
## [21] "Greplyng" "Nipdraugmose"
## [23] "Lavhei" "Sylmose"
## [25] "Torntvebladmose" "Fjellvåk"
## [27] "Jerv" "Kongeørn"
## [29] "Jaktfalk" "Fjellrev"
Dette er hovedsakelig en liste over arter eller artsgrupper. Untaket er lavhei og kanskje vier alpint belte. Det er litt spesielt at tamrein er med i denne indeksen/indikatoren siden vi har gjort en separat vurdering om at metoden bak denne indikatoren ikke er en som vi vil bruke selv.
Spesifiser hele landarealet til Norge, samt de tre regionene, som NIunits:
myNIunits <- c(allArea = T, parts = T, counties = F)
Inkludrer kun BSunits (kommuner) som har mer en 205 fjell:
myPartOfTotal <- .2
Siden denne opperasjonen tar litt tid så lagrer jeg outputen på server og henter det tilbake etterpå, så slipper jeg å kjøre gjennom hver gang.
fjell_assemeble <- NIcalc::assembleNiObject(
inputData = fjelldat,
predefNIunits = myNIunits,
partOfTotal = myPartOfTotal,
indexType = "ecosystem",
part = "ecosystem",
total = "terrestrial")
saveRDS(fjell_assemeble, "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellAssemble.rds")
rm(fjell_assemeble)
Henter tilbake datasettet
#fjell_assemeble <- readRDS("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellAssemble.rds")
fjell_assemeble <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellAssemble.rds")
Reduserer antall funsksjonelle grupper i datasettet for å ungå at noen indikatorer får veldig stor vekt pga at det er få indikatorer i den aktuelle gruppa. Mer spesifikt så slår i her sammen generalister og spesialister på hvert trofiske nivå.
xxx <- yyy <- fjell_assemeble$indicators$functionalGroup
xxxId <- yyyId <- fjell_assemeble$indicators$functionalGroupId
yyy[xxxId %in% c(1,2)] <- "Mellompredator"
yyyId[xxxId %in% c(1,2)] <- 1
yyy[xxxId %in% c(6,7)] <- "Primærprodusent"
yyyId[xxxId %in% c(6,7)] <- 6
yyy[xxxId %in% c(8,9)] <- "Topp-predator"
yyyId[xxxId %in% c(8,9)] <- 8
fjell_assemeble$indicators$functionalGroup <- yyy
fjell_assemeble$indicators$functionalGroupId <- yyyId
mountainDiagnostics <- NIcalc::imputeDiagnostics(x = fjell_assemeble,
nSim = 10,
transConst = 0.01,
maxit = 20)
mountainDiagnostics$diagnostics$mdDistribuiton
Det er noen datahull hvor standardsmetoden har imputerert data. (For 2019 er det kun snakk om vier og lavhei, så det kan være bedre for ØT-analysene å bare ekskludere disse.) Denne operasjonen tar en del tid, så jeg lagrer outputen på server og henter den tilbake (så slipper jeg å kjøre funksjonen på nytt hver gang).
fjell_imputes <- NIcalc::imputeData(x = fjell_assemeble,
nSim = 1000,
transConst = 0.01,
maxit = 20,
printFlag = TRUE)
saveRDS(fjell_imputes, "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-imputes.rds")
Henter tilbake datasettet
#fjell_imputes <- readRDS("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-imputes.rds")
fjell_imputes <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-imputes.rds")
mountainImputesCorr <- NIcalc::impStand(x = fjell_assemeble,
imputations = fjell_imputes)
Her regner jeg ut fjellindeksen tilsvarende den ‘offisielle’ versjonen
fjellIndex <- NIcalc::calculateIndex(
x = fjell_assemeble,
imputations = mountainImputesCorr,
nsim = 1000,
awBSunit = "Fjell",
fids = F, # should fidelities be ignored in
# the calculation of Wi?
tgroups = TRUE, # should grouping of indicators
# into trophic and key indicator
# groups be ignored
keys = "specialWeight", #"ignore",
w = 0.5
)
summary(fjellIndex$wholeArea)
## 2.5% median 97.5% displacement
## wholeArea 1990 0.5836179 0.6254757 0.6602899 -0.03283992
## wholeArea 2000 0.5534118 0.5921241 0.6242417 -0.02849441
## wholeArea 2010 0.5611825 0.6016798 0.6364247 -0.01743439
## wholeArea 2014 0.5154922 0.5734827 0.6127839 -0.03733600
## wholeArea 2019 0.5042129 0.5623335 0.5942140 -0.03486288
Disse tallene stemmer med det som presenteres på nettsiden. Tallene blit litt forskjellige hver gang pga permuteringen.
plot(fjellIndex$wholeArea)
La oss se hvordan denne standardmetoden har vektlagt de ulike indikatorene
wgt <- NIcalc::summaryWeights(fjellIndex$wholeArea)
wgt <- wgt[,ncol(wgt)]
wgt <- sort(wgt, decreasing = T)
par(mar=c(12,5,0.5,0))
barplot(wgt, las=2)
Fjellindeksen er i praksis oppgjort av kun rundt 6 indikatorer. De fire første er såkalt nøkkelindikatorer som har blitt gitt ekstra vekt. Rypeartene er ikke nøkkelarter, men får en del vekt, trolig fordi dataene har så god romlig dekning. Merk at 5 av de 6 høyest vektede indikatorene er med i ØT-fjell som egne, separete indikatorer fra før.
wgtT <- NIcalc::summaryWeights(fjellIndex$wholeArea, group = "troph")
wgtT<- wgtT[,ncol(wgtT)]
wgtT <- sort(wgtT, decreasing = T)
par(mar=c(12,5,0.5,0))
barplot(wgtT, las=2)
Her ser vi igjen hvor viktige nøkkelartene er i utregningen av indeksen. Standardmetoden vekter også ned indikatorer som representerer trofiske grupper hvor det allerede finnes flere andre indikatorer (lik vekt til alle trofiske nivå). Indeksen er med andre ord designet for å gi et bilde av tilstanden på et gitt areal, og ikke tilstanden til artene som inngår som indikatorer.
Standardmetoden er slik at skaleringen foregår for hver ICunit basert på indikator verdien og referanseverdien for dette arealet. Vi ønsker kanskje egentlig at skaleringen skal foregå bare for hver region (NIunits: e,s,w,n,c), samt at for totalen så bør den bare være felles for hele landet.
Noen indikatorer har bare en referanseverdi som er felles for hele Norge, slik som blåstrupe
DT::datatable(fjelldat$referenceValues$referenceValues[fjelldat$referenceValues$referenceValues$indName=="Blåstrupe", c("indName", "ICunitName", "expectedValue")])
Andre indikatorer har ulike referanseverdier for ulike arealenheter. Eksempelvis fjellrev:
DT::datatable(fjelldat$referenceValues$referenceValues[fjelldat$referenceValues$referenceValues$indName=="Fjellrev", c("indName", "ICunitName", "expectedValue")])
Vi har tenkt at at indikatorene og indeksen bør være indikatorfokusert og ikke areal-fokusert (slik NI er i dag). Vi så for oss at man kunne ta å summere (eller ta gjennomsnitt) av referanseverdien i alle ICunits i samme region og dele den summerte indikatorverdien på denne. For den skalerte indikatorverdien for hele landet kunne man summere referanseverdien i absolutt alle ICunits og dele den summerte indikatorverdien på denne.
Dette er ikke en funskjon i NIcalc og det krever en del utviklingsarbeid uten at det er opplagt at det i det hele tat er mulig. Det er for eksempel ikke opplagt hvordan man kan summer opp indikatorverdier og referanseverdier som er tettheter og ikke rene bestandstall. Vi beholder derfor den opprinnelige metoden for å bestemme hvor skaleringen skal foregå (dvs på ICunit-nivå).
Først vil jeg fjerne to indikatorer fra indikatorsettet - vier og lavhei. Disse har ikke data fra 2019 og siden vi ikke hovedsakelig er opptatt av tidserien og trenden i indeksen over tid, så seg jeg ingen grunn til å gjøre imputeringer. Jeg tror det letteste er å hente dataene på nytt fra databasen.
Her er listen over indikatorene jeg ønsker å ha med
(inds_red <- inds[!inds %in% c("Vier alpint belte", "Lavhei")])
## [1] "Heilo" "Boltit"
## [3] "Snøspurv" "Steinskvett"
## [5] "Heipiplerke" "Fjellerke"
## [7] "Ravn" "Blåstrupe"
## [9] "Ringtrost" "Lappspurv"
## [11] "Fjellrype" "Villrein"
## [13] "Lirype" "smågnagere - fjellbestander"
## [15] "Tamrein" "Fjellfiltmose"
## [17] "Issoleie" "Praktdraugmose"
## [19] "Fjellvalmue" "Greplyng"
## [21] "Nipdraugmose" "Sylmose"
## [23] "Torntvebladmose" "Fjellvåk"
## [25] "Jerv" "Kongeørn"
## [27] "Jaktfalk" "Fjellrev"
Her tar jeg inn bare de indikatorene som har data fra 2019 og så bryr jeg meg heller ikke om å lage en tidsserie, så jeg tar bare dataene fra 2019
fjelldat_red <- NIcalc::importDatasetApi(
username = myUser,
password = myPwd,
indic = inds_red,
year = c("1990","2000","2010","2014","2019"))
saveRDS(fjelldat_red, "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell_red_1990-2019.rds")
rm(fjelldat_red)
rm(myPwd)
Laster tilbake datasettet:
#fjelldat_red <- readRDS("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell_red_1990-2019.rds")
fjelldat_red <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell_red_1990-2019.rds")
Her må jeg velge indexType thematic siden jeg bruker mitt eget utvalg av indikatorer. Hvis ikke blir output en liste og ikke en niInput. Videre så skrur jeg av denne ekskluderingen av kommuner (BSunits) som har mindre en 20% fjellareal, siden jeg uansett vekter de skalerte indikatorene basert på total fjellareal senere.
fjell_assemeble_red <- NIcalc::assembleNiObject(
inputData = fjelldat_red,
predefNIunits = myNIunits,
indexType = "thematic"
)
saveRDS(fjell_assemeble_red, "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellAssemble_red_v2.rds")
rm(fjell_assemeble_red)
#fjell_assemeble_red <- readRDS("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellAssemble_red_v2.rds")
fjell_assemeble_red <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellAssemble_red_v2.rds")
Dette trenger jeg ikke gjøre nå ettersom jeg ikke skal bruke de trofiske nivåene til å bestemme vekting.
Det er ingen manglende data for 2019, men jeg kan gjøre imputering alikevel siden det er noen hull i de andre årene
mountainDiagnostics <- NIcalc::imputeDiagnostics(x = fjell_assemeble,
nSim = 10,
transConst = 0.01,
maxit = 20)
mountainDiagnostics_v2$diagnostics$mdDistribuiton
Det er noen datahull hvor standardsmetoden har imputerert data. Denne operasjonen tar en del tid, så jeg lagrer outputen på server og henter den tilbake (så slipper jeg å kjøre funksjonen på nytt hver gang).
fjell_imputes_v2 <- NIcalc::imputeData(x = fjell_assemeble_red,
nSim = 100, # litt lite med 100, men det er nok i dette tilfellet tror jeg
transConst = 0.01,
maxit = 20,
printFlag = TRUE)
saveRDS(fjell_imputes_v2, "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-imputes_v2.rds")
Henter tilbake datasettet
#fjell_imputes_v2 <- readRDS("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-imputes_v2.rds")
fjell_imputes_v2 <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjell-imputes_v2.rds")
mountainImputesCorr_v2 <- NIcalc::impStand(x = fjell_assemeble_red,
imputations = fjell_imputes_v2)
fjellIndex2 <- NIcalc::calculateIndex(
x = fjell_assemeble_red,
nsim = 10000,
awBSunit = "Fjell",
fids = F, # should fidelities be ignored in
# the calculation of Wi?
tgroups = F, # should grouping of indicators
# into trophic and key indicator
# groups be ignored
keys = "ignore",
awbs=TRUE # arealvekting basert på fjellareal i hver kommune
)
saveRDS(fjellIndex2, "P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellindeks2.rds")
fjellIndex2 <- readRDS("P:/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/fjell2021/data/naturindeks/NIfjellindeks2.rds")
Så lager jeg er dataframe med resultatene og eksporterer csv.
NIunits <- names(fjellIndex2)
myYears <- c("1990","2000","2010","2014","2019")
IND_samp <- data.frame(reg = rep(NIunits, each=length(myYears)*10000),
year = rep(myYears, each=10000, times= length(NIunits)),
val = NA)
for(i in 1:length(NIunits)){
for(n in 1:length(myYears)){
#print(i)
#print(n)
temp <- fjellIndex2[[i]][[n]]$index
IND_samp$val[IND_samp$reg == NIunits[i] &
IND_samp$year == myYears[n]] <- temp
}
}
IND_samp$reg[IND_samp$reg=="wholeArea"] <- "Norge"
head(IND_samp)
write.csv(IND_samp, "../output/indicator_values/fjellindeks.csv", row.names = F)
Her er verdiene for hele landet og for regionene.
tbl <- as.data.frame(summary(fjellIndex2$wholeArea))
tbl2 <- as.data.frame(summary(fjellIndex2$E))
tbl3 <- as.data.frame(summary(fjellIndex2$S))
tbl4 <- as.data.frame(summary(fjellIndex2$W))
tbl5 <- as.data.frame(summary(fjellIndex2$C))
tbl6 <- as.data.frame(summary(fjellIndex2$N))
NIfjellTbl <- rbind(tbl,
tbl2,
tbl3,
tbl4,
tbl5,
tbl6)
NIfjellTbl <- NIfjellTbl[,-4]
NIfjellTbl$reg = row.names(NIfjellTbl)
row.names(NIfjellTbl) <- seq(1,
nrow(NIfjellTbl),
1)
NIfjellTbl <- tidyr::separate(data=NIfjellTbl,
col = "reg",
into = c("reg",
"year"),
extra = "warn",
sep = " ",
remove=F)
NIfjellTbl$reg[NIfjellTbl$reg=="wholeArea"] <- "Norge"
names(NIfjellTbl)[names(NIfjellTbl)=="2.5%"] <- "low"
names(NIfjellTbl)[names(NIfjellTbl)=="median"] <- "med"
names(NIfjellTbl)[names(NIfjellTbl)=="97.5%"] <- "upp"
NIfjellTbl
## low med upp reg year
## 1 0.6568163 0.6906784 0.7232201 Norge 1990
## 2 0.6493923 0.6822344 0.7104447 Norge 2000
## 3 0.5945051 0.6413823 0.6856400 Norge 2010
## 4 0.6165059 0.6517108 0.6867010 Norge 2014
## 5 0.6183913 0.6524971 0.6864854 Norge 2019
## 6 0.6619410 0.7109596 0.7528823 E 1990
## 7 0.6157544 0.6668465 0.7115403 E 2000
## 8 0.5782963 0.6467766 0.7080557 E 2010
## 9 0.5782909 0.6274869 0.6774450 E 2014
## 10 0.5814563 0.6298341 0.6772609 E 2019
## 11 0.6790127 0.7310074 0.7789768 S 1990
## 12 0.5797079 0.6339178 0.6795633 S 2000
## 13 0.5391887 0.6122203 0.6778486 S 2010
## 14 0.5449940 0.5987425 0.6522464 S 2014
## 15 0.5557150 0.6073467 0.6594184 S 2019
## 16 0.6551510 0.7082790 0.7556510 W 1990
## 17 0.6204444 0.6732617 0.7204185 W 2000
## 18 0.5807376 0.6432760 0.7012042 W 2010
## 19 0.5858408 0.6358338 0.6847991 W 2014
## 20 0.5838044 0.6348691 0.6817723 W 2019
## 21 0.6661374 0.7146082 0.7586932 C 1990
## 22 0.6741869 0.7282261 0.7717091 C 2000
## 23 0.6258622 0.6905549 0.7488423 C 2010
## 24 0.6231656 0.6744452 0.7252280 C 2014
## 25 0.6249088 0.6747764 0.7232771 C 2019
## 26 0.6116576 0.6679186 0.7189445 N 1990
## 27 0.6344015 0.6829477 0.7218567 N 2000
## 28 0.5593559 0.6288349 0.6911869 N 2010
## 29 0.6131231 0.6638301 0.7128359 N 2014
## 30 0.6162245 0.6642395 0.7106803 N 2019
Indeksverdien blir en del høyere med denne metoden: 0.65 istedet for 0.56. (Kommentar: jeg fikk samme verdi 0.65 om jeg forstsatt hadde med lavhei og vier-indikatorene, så det hadde liten betydning i realiteten)
La oss se hvordan denne custom-metoden har vektlagt de ulike trofiske nivåene, derette de enklete indikatorene:
(wgtT_red <- NIcalc::summaryWeights(fjellIndex2$wholeArea, group = "troph"))
## Våtmark Skog wholeArea 1990 wholeArea 2000
## Mellompredator generalist 0 0 0.3436408 0.37969081
## Primærprodusent spesialist 0 0 0.1697326 0.16039454
## Topp-predator spesialist 0 0 0.1023101 0.09678519
## Plantespisere og filterspisere 0 30 0.1120618 0.10585468
## Topp-predator generalist 0 45 0.1697977 0.16031884
## Key indicator 80 25 0.1024570 0.09695594
## wholeArea 2010
## Mellompredator generalist 0.44071321
## Primærprodusent spesialist 0.14463600
## Topp-predator spesialist 0.08737414
## Plantespisere og filterspisere 0.09535333
## Topp-predator generalist 0.14433144
## Key indicator 0.08759188
Her skjer det noe rart. Jeg få bare ut deler av tabellen. Funksjonen fungere kanskje bare for type=ecosystem. Jeg kan regne ut vektene selv.
wgt2 <- NIcalc::summaryWeights(fjellIndex2$wholeArea, indicData=T)
DT::datatable(
wgt2,
extensions = "FixedColumns",
options = list(
scrollX = TRUE,
scrollY=T,
pageLength = 10
))
Fra datasettet over så bare summerer jeg vektene per funskjonelle gruppe, samt teller opp antall indikatorer per gruppe.
par(mar=c(12,5,1,0),
mfrow=c(1,2))
barplot(
tapply(
wgt2$`wholeArea 2019`,
wgt2$functionalGroup,
FUN = sum),
las=2,
ylab="Akkumulert indikatorvekt"
)
barplot(
tapply(
wgt2$`wholeArea 2019`,
wgt2$functionalGroup,
FUN = length),
las=2,
ylab="Antall indikatorer"
)
wgt2.1 <- NIcalc::summaryWeights(fjellIndex2$wholeArea, indicData=F)
wgt2.1 <- wgt2.1[,ncol(wgt2.1)]
wgt2.1 <- sort(wgt2.1, decreasing = T)
par(mar=c(12,5,0.5,0))
barplot(wgt2.1, las=2)
Nå var det mye mer jevnt. Ulikhetene nå skyldes to ting. 1) Arealvekting basert på fjellareal i hver kommune. Denne kan det også være lurt å beholde tror jeg. Se under oversikten ‘Index uten awbs’ under for hvordan indeksen ser ut da. 2) Ulik dekningsgrad - ikke alle indikatorene finnes eller har data fra alle BSunits (kommuner). Denne vektingen kan det nok være lurt å beholde. Merk at vektingen blir ulik for utregningen av de regionale indeksene. Med andre ord, tamrein har større vekt i nord enn i resten av landet (se under).
reinWgtW <- NIcalc::summaryWeights(fjellIndex2$W, indicData=T)
reinWgtN <- NIcalc::summaryWeights(fjellIndex2$N, indicData=T)
reinWgtW <- subset(reinWgtW, rownames(reinWgtW) %in% "Tamrein")
reinWgtN <- subset(reinWgtN, rownames(reinWgtN) %in% "Tamrein")
rein <- c(reinWgtW$`W 2019`, reinWgtN$`N 2019`)
names(rein) <- c("Tamrein Vest", "Tamrein Nord")
barplot(rein, ylab = "Indikatorvekt")
source("indicator_plots.R")
eval(parse("indicator_plots.R", encoding="UTF-8"))
fjellIndex3 <- NIcalc::calculateIndex(
x = fjell_assemeble,
imputations = mountainImputesCorr,
nsim = 500,
awBSunit = "Fjell",
fids = F, # should fidelities be ignored in
# the calculation of Wi?
tgroups = F, # should grouping of indicators
# into trophic and key indicator
# groups be ignored
keys = "ignore",
awbs=FALSE
)
La oss se hvordan denne custom-metoden har vektlagt de ulike indikatorene
wgt3 <- NIcalc::summaryWeights(fjellIndex3$wholeArea)
wgt3 <- wgt3[,ncol(wgt3)]
wgt3 <- sort(wgt3, decreasing = T)
par(mar=c(12,5,0.5,0))
barplot(wgt3, las=2)