1
$\begingroup$

As my subject suggests, I am running GAMs (in package mgcv) with sighting data as dependent and a variety of environmental and time parameters as predictor variables. I am running a hurdle model, so am modeling the presence absence data first with a quasibionomial distribution and then positive abundance model with gamma. I am selecting the best model using a forward step-wise approach. Each potential variable is added individually and the term resulting in the lowest GCV and highest deviance explained is then included in the next step. The trouble arises when the added variable in the model with the lowest GCV and highest deviance explained isn't actually significant.

Data:

structure(list(Grid_ID = structure(c(38L, 39L, 53L, 54L, 71L, 
    72L), .Label = c("1", "1,000", "1,001", "1,008", "1,009", "1,010", 
    "1,011", "1,012", "1,013", "1,014", "1,015", "1,016", "1,022", 
    "1,023", "1,024", "1,025", "1,026", "1,027", "1,028", "1,029", 
    "1,034", "1,035", "1,036", "1,037", "1,039", "1,040", "1,045", 
    "1,046", "1,047", "1,048", "1,053", "1,054", "1,055", "10", "100", 
    "101", "103", "104", "105", "106", "107", "108", "109", "11", 
    "110", "118", "119", "12", "121", "122", "123", "125", "126", 
    "127", "128", "129", "13", "130", "131", "132", "133", "14", 
    "141", "142", "143", "147", "148", "15", "150", "151", "152", 
    "153", "154", "155", "156", "157", "158", "159", "160", "161", 
    "162", "163", "167", "168", "169", "172", "173", "174", "175", 
    "176", "177", "178", "179", "180", "181", "182", "183", "184", 
    "185", "188", "189", "190", "194", "195", "196", "197", "198", 
    "199", "2", "20", "200", "201", "202", "203", "204", "205", "206", 
    "207", "209", "21", "210", "211", "217", "218", "219", "22", 
    "220", "221", "222", "223", "224", "225", "226", "227", "228", 
    "229", "23", "230", "231", "233", "234", "235", "236", "237", 
    "239", "246", "247", "248", "249", "25", "250", "252", "253", 
    "254", "255", "256", "257", "258", "259", "26", "260", "261", 
    "262", "267", "268", "269", "27", "270", "271", "272", "273", 
    "274", "275", "276", "277", "278", "279", "28", "280", "281", 
    "282", "285", "286", "287", "288", "289", "29", "290", "291", 
    "292", "293", "294", "295", "296", "297", "298", "299", "3", 
    "300", "301", "302", "303", "305", "306", "307", "308", "309", 
    "310", "311", "312", "313", "314", "315", "316", "317", "318", 
    "319", "320", "321", "326", "327", "328", "329", "330", "331", 
    "332", "333", "334", "335", "336", "337", "339", "340", "341", 
    "343", "344", "345", "346", "347", "348", "349", "350", "351", 
    "352", "355", "356", "357", "36", "360", "361", "362", "363", 
    "364", "365", "366", "367", "368", "369", "37", "372", "373", 
    "374", "376", "377", "378", "38", "380", "381", "382", "383", 
    "384", "385", "386", "39", "391", "392", "396", "397", "398", 
    "399", "4", "40", "400", "401", "402", "408", "409", "41", "410", 
    "412", "413", "414", "415", "416", "417", "42", "423", "424", 
    "425", "43", "430", "431", "432", "433", "434", "44", "441", 
    "442", "443", "444", "446", "447", "448", "449", "45", "450", 
    "451", "458", "459", "46", "460", "461", "462", "463", "466", 
    "467", "470", "471", "472", "473", "474", "475", "476", "484", 
    "485", "486", "487", "488", "489", "490", "491", "492", "495", 
    "496", "497", "498", "499", "5", "500", "501", "513", "514", 
    "515", "516", "517", "518", "523", "524", "525", "526", "527", 
    "528", "529", "54", "541", "542", "543", "544", "545", "55", 
    "550", "551", "552", "553", "554", "56", "569", "57", "570", 
    "571", "572", "573", "577", "578", "579", "58", "580", "581", 
    "582", "599", "60", "600", "601", "602", "603", "604", "605", 
    "606", "607", "608", "609", "61", "610", "62", "626", "627", 
    "628", "629", "63", "631", "632", "633", "634", "635", "636", 
    "637", "638", "639", "64", "653", "654", "655", "656", "657", 
    "658", "659", "661", "662", "663", "664", "665", "666", "667", 
    "668", "669", "670", "671", "672", "673", "687", "688", "689", 
    "690", "691", "692", "696", "697", "698", "699", "7", "700", 
    "701", "702", "703", "704", "705", "716", "717", "718", "719", 
    "720", "721", "722", "723", "724", "725", "726", "727", "728", 
    "739", "74", "740", "741", "745", "746", "747", "748", "749", 
    "75", "750", "751", "752", "753", "754", "764", "765", "766", 
    "767", "768", "769", "77", "770", "771", "772", "773", "78", 
    "782", "783", "784", "788", "789", "79", "790", "798", "799", 
    "8", "80", "800", "801", "804", "805", "81", "812", "813", "814", 
    "815", "816", "819", "82", "820", "821", "827", "828", "829", 
    "83", "830", "831", "833", "834", "835", "836", "84", "842", 
    "843", "844", "845", "846", "849", "85", "850", "851", "852", 
    "853", "854", "860", "861", "862", "863", "864", "869", "870", 
    "871", "872", "873", "874", "88", "881", "882", "883", "884", 
    "885", "886", "89", "890", "891", "892", "893", "894", "9", "902", 
    "903", "904", "905", "906", "908", "909", "910", "911", "912", 
    "922", "923", "924", "925", "926", "927", "928", "929", "930", 
    "940", "941", "942", "943", "944", "945", "946", "947", "948", 
    "957", "958", "959", "96", "960", "961", "962", "963", "964", 
    "965", "966", "97", "976", "977", "978", "979", "980", "981", 
    "982", "983", "984", "992", "993", "994", "995", "996", "997", 
    "998", "999"), class = "factor"), Lat = c(56.85614521, 56.85582097, 
    56.90062505, 56.90024495, 56.94504641, 56.94461032), Long = c(-153.5604891, 
    -153.4783612, -153.4777153, -153.3954873, -153.3947378, -153.3124098
    ), Er_Pres = c(0L, 0L, 0L, 0L, 0L, 0L), Er_Count = c(0L, 0L, 
    0L, 0L, 0L, 0L), Mn_Pres = c(0L, 0L, 0L, 0L, 0L, 0L), Mn_Count = c(0L, 
    0L, 0L, 0L, 0L, 0L), Bp_Pres = c(0L, 0L, 0L, 0L, 0L, 0L), Bp_Groups = c(0L, 
    0L, 0L, 0L, 0L, 0L), Bp_Count = c(0L, 0L, 0L, 0L, 0L, 0L), Mn_Groups = c(0L, 
    0L, 0L, 0L, 0L, 0L), Month = c(8L, 8L, 8L, 8L, 8L, 8L), Year = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L), .Label = c("1997", "1998", "1999", "2000", 
    "2001", "2002", "2003", "2004", "2005", "2006", "2007", "2008", 
    "2009", "2010", "2011", "2012", "2013"), class = "factor"), Er_Count_Density = c(0, 
    0, 0, 0, 0, 0), Mn_Count_Density = c(0L, 0L, 0L, 0L, 0L, 0L), 
        Bp_Count_Density = c(0, 0, 0, 0, 0, 0), Bathymetry = c(-64.1858139, 
        -76.11354065, -92.14147949, -90.60312653, -80.86385345, -71.55316162
        ), Grid_Area = c(25, 25, 25, 25, 25, 25), DFS = c(3.688910639, 
        6.807817092, 4.233185446, 9.199096676, 7.19354409, 5.153224038
        ), Slope = c(0.15274878, 0.13670446, 0.38316911, 0.08646853, 
        0.1305674, 0.20038579), DOY = c(244L, 244L, 244L, 244L, 244L, 
        244L), Quarter = structure(c(3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1", 
        "2", "3", "4"), class = "factor"), SST = c(12.6000004, 12.75, 
        12.8999996, 13.0500002, 12.5249996, 12.6000004), chla = c(0.53747, 
        0.53747, 0.53747, 0.534891, 0.58174, 0.696316), Depth = c(64.1858139, 
        76.11354065, 92.14147949, 90.60312653, 80.86385345, 71.55316162
        )), .Names = c("Grid_ID", "Lat", "Long", "Er_Pres", "Er_Count", 
    "Mn_Pres", "Mn_Count", "Bp_Pres", "Bp_Groups", "Bp_Count", "Mn_Groups", 
    "Month", "Year", "Er_Count_Density", "Mn_Count_Density", "Bp_Count_Density", 
    "Bathymetry", "Grid_Area", "DFS", "Slope", "DOY", "Quarter", 
    "SST", "chla", "Depth"), row.names = c(NA, 6L), class = "data.frame")

GAM Model

gam2g<-gam(Er_Pres~ te(Lat, Long,k=10,bs='tp')+ Year,family=quasibinomial(link='logit'),data=Data,gamma=1.4,offset(Grid_Area))

Results from summary(gam2g)

Family: quasibinomial 
Link function: logit 

Formula:
Er_Pres ~ te(Lat, Long, k = 10, bs = "tp") + Year

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)   -1817.84 7785870.42       0        1
Year1998        -14.06 9210994.35       0        1
Year1999        462.72 7785870.12       0        1
Year2000        462.62 7785870.12       0        1
Year2001        463.56 7785870.12       0        1
Year2002        463.58 7785870.12       0        1
Year2003        463.18 7785870.12       0        1
Year2004        463.17 7785870.12       0        1
Year2005        463.12 7785870.12       0        1
Year2006         98.36 8018436.07       0        1
Year2007         58.53 8229756.26       0        1
Year2008        464.76 7785870.12       0        1
Year2009        465.33 7785870.12       0        1
Year2010        464.53 7785870.12       0        1
Year2011        463.94 7785870.12       0        1
Year2012        463.11 7785870.12       0        1
Year2013        464.09 7785870.12       0        1

Approximate significance of smooth terms:
               edf Ref.df     F p-value    
te(Lat,Long) 49.79  51.97 13.79  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.193   Deviance explained = 43.5%
GCV = 2.3191  Scale est. = 7.5651    n = 12458

This model produces the lowest GCV with the highest deviance explained versus all other model options at this step, but none of the year coefficients are significant. Am I justified in selecting year at this step? The results seem meaningless.

$\endgroup$
  • $\begingroup$ looking at the year parameter estimates, they're all about the same except for 1997, 1998, 2006 and 2007. You could estimate coefficients for just those years (recoding Data$strange_years <- factor(ifelse(as.character(Data$Year) %in% c("1997", "1998", "2006", "2007"), as.character(Data$Year), "normal year"))). How many observations do you have per year? $\endgroup$ – Gregor Oct 28 '14 at 19:20
  • $\begingroup$ Why do you want year as a factor? The huge standard errors suggest real problems with this model as it stands. You could try to convert Year to an integer and then use s(Year) in the model for example. If you do this, you can use select = TRUE, method = "REML" to turn on additional penalties that can shrink terms out of the model - which performs a sort of feature selection for you in a much more principled manner than stepwise selection. $\endgroup$ – Gavin Simpson Oct 28 '14 at 19:39
  • $\begingroup$ I'm including year as a factor because I want to see if there are temporal changes in the observations. I'm also exploring seasonal trends by looking at Quarter as a factor. $\endgroup$ – akbreezo Oct 28 '14 at 21:21
  • $\begingroup$ There is definitely less effort in the "strange years." There is only one survey in 1997 and 1998 and 3 each in 2006 and 2007. Every other year has at least 5. I thought I accounted for this variance in effort by including effort as an offset. Perhaps there is a better way to handle this? $\endgroup$ – akbreezo Oct 28 '14 at 21:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.