The Bowman Lab

Hello! My name is Melissa Hopkins. I just finished my first quarter as an undergraduate researcher in the Bowman Lab. The project I am working on involves the diversity of halophiles from the South Bay Saltworks lakes in San Diego. The South Bay Saltworks is an active solar salt harvesting facility that is part of the San Diego Bay National Wildlife Refuge. The goal of this project is to use 16S/18S community structure to identify microbial taxa that are currently poorly represented in existing genomes. We want to contribute to halophile genomes to learn more about halophiles, and this helps us decide which new genomes to sequence.

Aharan Oren’s open access review Microbial life at high concentrations” phylogenetic and metabolic diversity” explains the different classes and orders that contain halophiles, as well as the similarities of strategies used by these halophiles. Here, he uses the definition of halophiles as microbes that are able to tolerate 100 g/L salt but grow optimally at 50 g/L salt (seawater contains 35 g/L salt). Extreme halophiles (including the haloarchaea) are defined as growing best at salt concentrations of 2.5-5.2 M, and moderate halophiles as growing optimally at salt concentrations of 1.5-4.0 M. This assumes that the salt is sodium chloride however, different salts such as magnesium chloride can present additional challenges to life.

Halophiles are able to survive these high salt concentration environments in two different ways: pumping salt ions into their cells from the surrounding environment, or synthesizing organic solutes to match the concentration of their surrounding environment. Synthesizing organic solutes is more energetically expensive because it requires energy to make the high concentration of organic solutes needed, thus keeping salt ions out of the cytoplasm. But, this strategy is actually found widely across halophiles species. Pumping specific salt ions into their cells that don’t interfere with biological processes is less energetically expensive, but requires proteins in the cell to be specially adapted to high salt conditions. Because of this that strategy is not seen as often across the different species of halophiles. Different families and orders of halophiles use variations of these strategies to survive. There have been some new halophile species, such as Salinbacter, that are outliers, as in they use a different survival strategy then other halophiles they are related to. Sequencing many more halophile genomes will give us new information on how these adaptive strategies across different halophiles.

For this project, Jeff, Natalia and I spent the day sampling lakes from South Bay Saltworks on October 6. Out goal was to sample from 3 different points in each of several lakes of different salt concentrations, and to sample as many lakes as possible.

We started out at the lower salinity lakes to see how the equipment would function and get used to the sampling process. At each point, we took unfiltered samples using a peristaltic pump for respiration tests, bacterial abundance, measurements of photosynthetic efficiency and chlorophyll concentration, turbidity, and salinity.

We then placed a GF/F filter on the housing of the pump to collect a coarsely filtered sample for ion composition analysis, FDOM, and dissolved inorganic nutrients. Finally, we placed a 0.2 micron filter on the housing to collect bacteria, archaea, and phytoplankton for DNA analysis.

In all, we sampled 7 lakes: one low salt concentration lake, one medium salt concentration lake, 2 high salt concentration lakes, and 3 magnesium chloride lakes. Unfortunately, due to time constraints and the high viscosity of the magnesium chloride lakes (its like trying to filter maple syrup), we were only to sample from 1 point for each of the magnesium chloride lakes.

Natalia and I setting up for sampling on one of the lower salinity lakes.  You can see the large piles of harvested salt in the background.

One of the magnesium chloride lakes we sampled. Due to the high salt concentration and extremely high attraction between water and magnesium chloride these lakes have an oily texture and high viscosity, making them difficult to sample.

One of the saltier sodium chloride lakes that we sampled from.  The pink color comes from pigments in the different microorganisms that we’re studying.

Today we took the last sample of our Ecuador field effort, though we have a few days left in-country. Right now we are in the town of Mompiche, just down the coast from our second field site near Muisne. Tomorrow we’ll be sorting out gear and getting ready for a few days of meetings in Guayaquil. Then its time to fly home and start working up some data! I’m too tired to write a coherent post on the last two (intensive!) weeks of sampling, but here’s a photo summary of our work in the Cayapas-Mataje Ecological Reserve, where we spent the bulk of our time. Enjoy!

Pelicans congregating along a front in the Cayapas-Mataje estuary.

The town of Valdez, near the mouth of the Cayapas-Mataje estuary.  The Reserve is right on the border with Columbia, and up until a few years ago Valdez had a reputation as a trafficking hub.  Drug trafficking is still a problem throughout the Reserve, but with the conflict with FARC more or less over I understand that tension in the local communities has gone way down.  Valdez seems okay, and the people we met there were friendly.

Another view of Valdez.

Shrimp farm in the Cayapas-Mataje Ecological Reserve.  You can’t build a new farm in the Reserve, but old shrimp farms were grandfathered in when the Reserve was created.

The R/V Luz Yamira, at its home port of Tambillo.  Tambillo was a vibrant, friendly little town where we spent a bit of time.  The town is struggling to hold onto its subsistence fishing lifestyle in the face of declining fish stocks.

ADCP, easy as 1-2-3

Birds of a feather…

Morning commute from the city of San Lorenzo.

The Cayapas-Mataje Ecological Reserve has the tallest mangrove trees in the world.

I took this picture just for the good people at UCSD risk management.

Team Cayapas-Mataje.  From left; Jessie, Natalia, and Jeff.  We are standing in front of Jessie’s house in Tambillo.  Many thanks to Jessie and his wife for letting us stay a night and get away from the craziness of San Lorenzo!

A very full car ready to head to Muisne.  Its a good thing Natalia and I are both fairly short.

A large shrimp farm near Muisne.  The area around Muisne has been almost entirely deforested for shrimp aquaculture.  By comparing this area with the more pristine Cayapas-Mataje, we hope to better understand the biogeochemical consequences of coastal land-use change.

As any reader of this blog will know, most of the research in the Bowman Lab is focused on polar microbial ecology. Although focusing a research program on a set of geographically-linked environments does have advantages, primarily the ability to spend more time thinking in depth about them, there is I think something lost with this approach. Insights are often gained by bringing a fresh perspective to a new study area, or applying lessons learned in such an area to places that one has studied for years. With this in mind lab member Natalia Erazo and I are launching a new field effort in coastal Ecuador. Natalia is an Ecuadorean native, and gets credit for developing the idea and sorting out the considerable logistics for this initial effort. Our first trip is very much a scouting effort, but will carry out some sampling in the Cayapas-Mataje Ecological Reserve and near the town of Muisne. Depending on funding we hope to return during the rainy season in January-February for a more intensive effort.

Dense coastal forest in the Cayapas-Mataje Ecological Reserve.  Photo Natalia Erazo.

Our primary objective is to understand the role of mangrove forests in coastal biogeochemical cycling. Mangroves are salt-tolerant trees that grow in tropical and sub-tropical coastal areas around the world. They are known to provide a range of positive ecosystem functions; serving as fish habitat, stabilizing shorelines, and providing carbon and nutrient subsidies to the coastal marine environment. Globally mangroves are under threat. The population density of many tropical coastal areas is increasing, and that inevitably leads to land-use changes (such as deforestation) and a loss of economic services – the social and economic benefits of an ecosystem – as economic activity ramps up. The trick to long-term sustainability is to maintain ecosystem services during economic development, allowing standards of living to increase in the short term without a long-term economic loss resulting from ecological failure (such as the collapse of a fishery or catastrophic coastal flooding). This is not easily done, and requires a much better understanding of what functions exactly, specific at-risk components of an ecosystem provide than we often have.

One particular land-use threat to the mangrove ecosystem is shrimp aquaculture. Mangrove forests in Ecuador and in other parts of the world have been deforested on a massive scale to make room for shrimp aquaculture ponds. In addition to scaling back any ecosystem functions provided by the mangrove forest, the shrimp aquaculture facilities are a source of nutrients in the form of excrement and excess feed. On this trip we will try to locate estuaries more and less perturbed by aquaculture. By comparing nutrient and carbon concentrations, sediment load, and microbial community structure between these areas, we will gain a preliminary understanding of what happens to the coastal ecosystem when mangroves are removed and aquaculture facilities are installed in their place.

Our first stop on this search will be San Lorenzo, a small city in the Cayapas-Mataje Ecological Reserve near the border with Columbia. I’m extremely excited to visit the Reserve, which has the distinction of hosting the tallest mangrove trees anywhere on Earth. We may some meager internet access in San Lorenzo, so I’ll try to update this blog as I’m able. Because of the remote nature of some of our proposed field sites we’ll have a Garmin InReach satellite messenger with us. We plan to leave the device on during our field outings, you can track our location in real time on the map below! The Garmin map interface is a bit cludgey; you should ignore the other Scripps Institution of Oceanography users listed on the side panel as I can’t seem to make them disappear.

Davos, Switzerland, site of the POLAR2018 conference.  Image from https://www.inghams.co.uk

With colleagues Maria Corsaro, Eric Collins, Maria Tutino, Jody Deming, and Julie Dinasquet I’m convening a session on polar microbial ecology and evolution at the upcoming POLAR2018 conference in Davos, Switzerland.  Polar2018 is shaping up to be a unique and excellent conference; for the first time (I think) the major international Arctic science organization (IASC) is joining forces with the major international Antarctic science organization (SCAR) for a joint meeting.  Because many polar scientists specialize in one region, and thus have few opportunities to learn from the other, a joint meeting will be a great opportunity to share ideas.

Here’s the full-text description for our session.  If your work is related to microbial evolution, adaption, or ecological function in the polar regions I hope you’ll consider applying!

The eastern oyster.  Denitrification never tasted so good!  Photo from http://dnr.sc.gov.

I’m happy to be co-author on a study that was just published by Ann Arfken, a PhD student at the Virginia Institute for Marine Science (VIMS).  The study evaluated the composition of the microbial community associated with eastern oyster Crassostrea virginica to determine if oysters play a role in denitrification.  Denitrification is an ecologically significant anaerobic microbial metabolism.  In the absence of oxygen certain microbes use other oxidized compounds as electron acceptors.  Nitrate (NO₃^–) is a good alternate electron acceptor, and the process of reducing nitrate to nitrite (NO₂), and ultimately to elemental nitrogen (N₂), is called denitrification.  Unfortunately nitrate is needed by photosynthetic organisms like plants and phytoplankton, so the removal of nitrate can be detrimental to ecosystem health.  Oxygen is easily depleted in the guts of higher organisms by high rates of metabolic activity, creating a niche for denitrification and other anaerobic processes.

Predicted relative abundance (paprica) as a function of measured (qPCR) relative abundance of nosZI genes.  From Arfken et al. 2017.

To evaluate denitrification in C. virginica, Arfken et al. coupled actual measurements of denitrification in sediments and oysters with an analysis of microbial community structure in oyster shells and digestive glands.  We then used paprica with a customized database to predict the presence of denitrification genes, and validated the predictions with qPCR.

I was particularly happy to see that the qPCR results agreed well with the paprica predictions for the nosZ gene, which codes for the enzyme responsible for reducing nitrous oxide (N₂O) to N₂.  I believe this is the first example of qPCR being used to validate metabolic inference – which currently lacks a good method for validation.  Surprisingly however (at least to me), denitrification in C. virginica was largely associated with the oyster shell rather than the digestive gland.  We don’t really know why this is.  Arfken et al. suggests rapid colonization of the oyster shell by denitrifying bacteria that also produce antibiotic compounds to exclude predators, but further work is needed to demonstrate this!

Artist’s rendition of the Dawn spacecraft approaching Ceres.  Original image can be found here.

I’m excited to be hosting the fall meeting of NASA’s Outer Planets Assessment Group (OPAG) here at Scripps Institution of Oceanography in September.  For planetary scientists at UCSD, SDSU, USD, and other institutions in the greater San Diego area, if you’ve never attended the OPAG meeting here’s your chance!  The meeting will be September 6 and 7 at the Samual H. Scripps Auditorium.  More details can be found here.

What are assessment groups, and what is OPAG specifically?  The assessment groups are an excellent NASA innovation to encourage dialogue within the scientific community, and between the scientific community and NASA HQ.  There’s usually a little tense dialogue – in a good way – between these two ends of the scientific spectrum.  I often wish NSF had a similar open-format dialogue with its user community!  The form of the OPAG meeting is 15 or 30 minute presentations on a range of topics relevant to the community.  These are often mission updates, planning updates for future missions, and preliminary results from the analysis of mission data.  NASA has quite a few assessment groups, ranging from the Small Body Assessment Group (SBAG – probably the AG with the catchiest acronym) to the Mars Assessment Group (MPAG).  OPAG covers everything in the solar system further from the sun than Mars.  If that covers your favorite planetary body, come and check it out!

It’s traditional to have a public evening lecture with the OPAG meeting.  For the upcoming meeting the lecture will be given at 7 pm on September 6 at the Samuel Scripps Auditorium by my friend and colleague Britney Schmidt, a planetary scientist from Georgia Tech and an expert on Europa and on Antarctic ice sheets.  Why and how one can develop that dual expertise will certainly be made clear in her talk.  There is no cost to attend, but an RSVP is recommended.  You can find more details and RSVP here.

Very often environmental datasets contain far fewer observations than we would like, and far more variables that might influence these observations.  This is the case for one of my projects in the Palmer LTER: I have annual observations of an environmental indicator (n = 19) and far more than 19 environmental conditions (predictors, p) that might influence the indicator.  This is a broad dataset (n >> p), and the problem can’t be solved with classic multiple regression (for which the number of observations must exceed the number of predictors).  Enter the R package glmnet.  Glmnet is an implementation of lasso, ridge, and elastic-net regression.  There are a limited number of glmnet tutorials out there, including this one, but I couldn’t find one that really provided a practical start to end guide.  Here’s the method that I came up with for using glmnet to select variables for multiple regression.  I’m not an expert in glmnet, so don’t take this as a definitive guide.  Comments are welcome!

First, let’s get some data.  The data are various climate indices with a 13- month time lag as well as sea ice extent, open water extent, and fractional open water extent.  There are a couple of years missing from this dataset, and we need to remove one more year for reasons that aren’t important here.

env.params <- read.csv('https://www.polarmicrobes.org/extras/env_params.csv', header = T, row.names = 1)
env.params <- env.params[row.names(env.params) != '1994',]

To keep things simple we’ll create our response variable by hand. Nevermind what the response variable is for now, you can read our paper when it comes out 🙂

response <- c(0.012, 0.076, 0.074, 0.108, 0.113, 0.154, 0.136, 0.183, 0.210, 0.043, 0.082, 0.092, 0.310, 0.185, 0.436, 0.357, 0.472, 0.631, 0.502)

One thing to note at this point is that, because we have so few observations, we can’t really withhold any to cross-validate the glmnet regression. That’s a shame because it means that we can’t easily optimize the alpha parameter (which determines whether glmnet uses lasso, elastic net, or ridge) as was done here. Instead we’ll use alpha = 0.9. This is to make the analysis a little bit elastic-netty, but still interpretable.

In my real analysis I was using glmnet to identify predictors for lots of response variables.  It was therefor useful to write a function to generalize the several commands needed for glmnet.  Here’s the function:

library(glmnet)

## The function requires a matrix of possible predictors, a vector with the response variable,
## the GLM family used for the model (e.g. 'gaussian'), the alpha parameter, and "type.measure".
## See the documentation on cv.glmnet for options.  Probably you want "deviance". 

get.glmnet <- function(predictors, response, family, alpha, type.measure){
  glmnet.out <- glmnet(predictors, response, family = family, alpha = alpha)
  glmnet.cv <- cv.glmnet(predictors, response, family = family, alpha = alpha, type.measure = type.measure)
  
  ## Need to find the local minima for lambda.
  lambda.lim <- glmnet.cv$lambda[which.min(glmnet.cv$cvm)]
  
  ## Now identify the coefficients that correspond to the lambda minimum.
  temp.coefficients.i <- coefficients(glmnet.out, s = lambda.lim)@i + 1 # +1 to account for intercept
  
  ## And the parameter names...
  temp.coefficients.names <- coefficients(glmnet.out, s = lambda.lim)@Dimnames[[1]][temp.coefficients.i]
  temp.coefficients <- coefficients(glmnet.out, s = lambda.lim)@x
  
  ## Package for export.
  temp.coefficients <- rbind(temp.coefficients.names, temp.coefficients)
  return(temp.coefficients)
}

Phew!  Okay, let’s try to do something with this.  Note that glmnet requires a matrix as input, not a dataframe.

response.predictors <- get.glmnet(data.matrix(env.params), response, 'gaussian', alpha, 'deviance')

>t(response.predictors)
      temp.coefficients.names temp.coefficients      
 [1,] "(Intercept)"           "0.496410195525042"    
 [2,] "ao.Apr"                "0.009282064516813"    
 [3,] "ao.current"            "-0.0214919836174853"  
 [4,] "pdo.Mar"               "-0.0568728879266135"  
 [5,] "pdo.Aug"               "-0.00881845994191182" 
 [6,] "pdo.Dec"               "-0.0321738320415234"  
 [7,] "ow.Dec"                "1.92231198892721e-06" 
 [8,] "ow.frac.current"       "0.207945851122607"    
 [9,] "ice.Sep"               "-2.29621552264475e-06"

So glmnet thinks there are 9 predictors.  Possible, but I’m a little suspicious.  I’d like to see this in the more familiar GLM format so that I can wrap my head around the significance of the variables.  Let’s start by building a null model of all the predictors.

null.model <- lm(response ~ ao.Apr + ao.current + pdo.Mar + pdo.Aug + pdo.Dec + ow.Dec + ow.frac.current + ice.Sep, data = env.params)

> summary(null.model)

Call:
lm(formula = response ~ ao.Apr + ao.current + pdo.Mar + pdo.Aug + 
    pdo.Dec + ow.Dec + ow.frac.current + ice.Sep, data = env.params)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.14304 -0.03352 -0.01679  0.05553  0.13497 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)      2.052e+00  1.813e+00   1.132   0.2841  
ao.Apr           2.226e-02  3.702e-02   0.601   0.5610  
ao.current      -3.643e-02  1.782e-02  -2.044   0.0681 .
pdo.Mar         -7.200e-02  3.215e-02  -2.240   0.0490 *
pdo.Aug         -1.244e-02  2.948e-02  -0.422   0.6821  
pdo.Dec         -6.072e-02  3.664e-02  -1.657   0.1285  
ow.Dec           3.553e-06  1.749e-06   2.032   0.0696 .
ow.frac.current  1.187e-01  4.807e-01   0.247   0.8100  
ice.Sep         -1.023e-05  8.558e-06  -1.195   0.2596  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09043 on 10 degrees of freedom
Multiple R-squared:  0.8579,	Adjusted R-squared:  0.7443 
F-statistic: 7.549 on 8 and 10 DF,  p-value: 0.002224

Hmmm… this looks messy.  There’s only one predictor with a slope significantly different from 0.  A high amount of variance in the original data is accounted for, but it smells like overfitting to me.  Let’s QC this model a bit by 1) checking for autocorrelations, 2) using AIC and relative likelihood to further eliminate predictors, and 3) selecting the final model with ANOVA.

## Check for autocorrelations using variance inflation factors

library(car)

> vif(null.model)
         ao.Apr      ao.current         pdo.Mar         pdo.Aug         pdo.Dec          ow.Dec ow.frac.current 
       1.610203        1.464678        1.958998        2.854688        2.791749        1.772611        1.744653 
        ice.Sep 
       1.510852

Surprisingly, all the parameters have acceptable vif scores (vif < 5).  Let’s proceed with relative likelihood.

## Define a function to evaluate relative likelihood.

rl <- function(aicmin, aici){
  return(exp((aicmin-aici) / 2))
}

## Construct multiple possible models, by adding parameters by order of abs(t value) in summary(null.lm).

model.1 <- lm(response ~ pdo.Mar, data = env.params)
model.2 <- lm(response ~ pdo.Mar + ao.current, data = env.params)
model.3 <- lm(response ~ pdo.Mar + ao.current + ow.Dec, data = env.params)
model.4 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec, data = env.params)
model.5 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep, data = env.params)
model.6 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + ao.Apr, data = env.params)
model.7 <- lm(response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + ao.Apr + pdo.Aug, data = env.params)

## Collect AIC scores for models.

model.null.aic <- AIC(null.model)
model.1.aic <- AIC(model.1)
model.2.aic <- AIC(model.2)
model.3.aic <- AIC(model.3)
model.4.aic <- AIC(model.4)
model.5.aic <- AIC(model.5)
model.6.aic <- AIC(model.6)
model.7.aic <- AIC(model.7)

## Identify the model with the lowest AIC score.

> which.min(c(model.1.aic, model.2.aic, model.3.aic, model.4.aic, model.5.aic, model.6.aic, model.7.aic, model.null.aic))
[1] 5

So model.5 has the lowest AIC.  We need to check the relative likelihood of other models minimizing information loss. Models with values < 0.05 do not have a significant likelihood of minimizing information loss and can be discarded.

> rl(model.5.aic, model.1.aic)
[1] 8.501094e-05
> rl(model.5.aic, model.1.aic)
[1] 8.501094e-05
> rl(model.5.aic, model.2.aic)
[1] 0.00143064
> rl(model.5.aic, model.3.aic)
[1] 0.002415304
> rl(model.5.aic, model.4.aic)
[1] 0.7536875
> rl(model.5.aic, model.6.aic)
[1] 0.4747277
> rl(model.5.aic, model.7.aic)
[1] 0.209965
> rl(model.5.aic, model.null.aic)
[1] 0.08183346

Excellent, we can discard quite a few possible mode here; model.1, model.2, and model.3.  Last we’ll use ANOVA and the chi-squared test to see if any of the models are significantly different from model.4, the model with the fewest parameters.

> anova(model.4, model.5, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     13 0.086169  1  0.012457   0.1704

> anova(model.4, model.6, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + 
    ao.Apr
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     12 0.083887  2   0.01474   0.3485

> anova(model.4, model.7, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec + ice.Sep + 
    ao.Apr + pdo.Aug
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     11 0.082276  3   0.01635   0.5347

> anova(model.4, null.model, test = 'Chisq')
Analysis of Variance Table

Model 1: response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec
Model 2: response ~ ao.Apr + ao.current + pdo.Mar + pdo.Aug + pdo.Dec + 
    ow.Dec + ow.frac.current + ice.Sep
  Res.Df      RSS Df Sum of Sq Pr(>Chi)
1     14 0.098626                      
2     10 0.081778  4  0.016849   0.7247

Okay, there’s a whole lot going on there, they important thing is that none of the models are significantly different from the model with the fewest parameters. There are probably some small gains in the performance of those models, but at an increased risk of over fitting.  So model.4 is the winner, and we deem pdo.Mar, ao.current, ow.Dec, and pdo.Dec to be the best predictors of the response variable.  For those of you who are curious, that’s the March index for the Pacific Decadal Oscillation, the index for the Antarctic Oscillation when the cruise was taking place, within-pack ice open water extent in December (the month prior to the cruise), and the Pacific Decadal Oscillation index in December.  Let’s take a look at the model:

> summary(model.4)

Call:
lm(formula = response ~ pdo.Mar + ao.current + ow.Dec + pdo.Dec, 
    data = env.params)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.177243 -0.037756 -0.003996  0.059606  0.114155 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.586e-02  5.893e-02  -0.609 0.552563    
pdo.Mar     -9.526e-02  2.208e-02  -4.315 0.000713 ***
ao.current  -3.127e-02  1.590e-02  -1.967 0.069308 .  
ow.Dec       4.516e-06  1.452e-06   3.111 0.007663 ** 
pdo.Dec     -8.264e-02  2.172e-02  -3.804 0.001936 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08393 on 14 degrees of freedom
Multiple R-squared:  0.8287,	Adjusted R-squared:  0.7797 
F-statistic: 16.93 on 4 and 14 DF,  p-value: 2.947e-05

You’ll recall that the null model accounted for 74 % of the variance in the response variable, our final model accounts for 78 % and we’ve dropped several parameters.  Not bad!  Note that we never could have gotten to this point however, without the holistic search provided by glmnet.

I can still hear some grumbling about over fitting however, so let’s try to address that by bootstrapping.  We are crazy data-limited with only 19 observations, but let’s iteratively hold three observations back at random, build a model with the remaining 16 (using our identified parameters), and see how well each model predicts the three that were held back.  I created a function to do this so that I could repeat this exercise for lots of different response variables.

## The function requires as input the response vector, a vector of the predictors,
## the data frame of predictors, and the number of iterations you want to run.

bootstrap.model <- function(response.vector, predictors, data, n){
  
  predictor.df <- as.data.frame(data[,predictors])
  colnames(predictor.df) <- predictors
  model.predictions <- matrix(ncol = 3, nrow = 0)
  colnames(model.predictions) <- c('year', 'prediction', 'real')
  
  ## How many observations you need to withhold is determined by how many
  ## iterations you want to run.  For example, for 1000 iterations of
  ## unique random subsets of 19 observations, you need to withold 3.
  
  x <- 0
  boot <- 0
  
  for(y in (length(response.vector) - 1):1){
    while(x < n){
      boot <- boot + 1
      x <- y ** boot
      print(x)
    }
  }
  
  ## Now build n new models.
  
  for(i in 1:n){
    print(i)
    predict.i <- sample.int(length(response.vector), boot)
    train.i <- which(!1:length(response.vector) %in% predict.i)
    
    temp.model <- lm(response.vector ~ ., data = predictor.df, subset = train.i)
    
    new.data <- as.data.frame(predictor.df[predict.i,])
    colnames(new.data) <- predictors
    
    temp.predict <- predict(temp.model, newdata = new.data, type = 'response')
    temp.actual <- response.vector[predict.i]
    
    temp.out <- matrix(ncol = 3, nrow = boot)
    try({temp.out[,1] <- names(response.vector)[predict.i]}, silent = T)
    temp.out[,2] <- temp.predict
    temp.out[,3] <- temp.actual
    
    model.predictions <- rbind(model.predictions, temp.out)
    
  }
  
  mean.predicted <- as.data.frame(tapply(as.numeric(model.predictions[,2]), as.character(model.predictions[,3]), mean))
  sd.predicted <- as.data.frame(tapply(as.numeric(model.predictions[,2]), as.character(model.predictions[,3]), sd))
  
  ## Make an awesome plot.

  plot(mean.predicted[,1] ~ as.numeric(row.names(mean.predicted)),
     ylab = 'Predicted response',
     xlab = 'Observed response',
     pch = 19)
  
  for(r in row.names(mean.predicted)){
    lines(c(as.numeric(r), as.numeric(r)), c(mean.predicted[r,1], mean.predicted[r,1] + sd.predicted[r,1]))
    lines(c(as.numeric(r), as.numeric(r)), c(mean.predicted[r,1], mean.predicted[r,1] - sd.predicted[r,1]))
  }
  
  abline(0, 1, lty = 2)
  
  mean.lm <- lm(mean.predicted[,1] ~ as.numeric(row.names(mean.predicted)))
  abline(mean.lm)
  print(summary(mean.lm))
  
  ## Capture stuff that might be useful later in a list.
  
  list.out <- list()
  list.out$model.predictions <- model.predictions
  list.out$mean.lm <- mean.lm
  list.out$mean.predicted <- mean.predicted
  list.out$sd.predicted <- sd.predicted
  
  return(list.out)
  
}

## And execute the function...

demo.boostrap <- bootstrap.model(response, c('pdo.Mar', 'ao.current', 'ow.Dec', 'pdo.Dec'), env.params, 1000)

Hey, that looks pretty good!  The dotted line is 1:1, while the solid line gives the slope of the regression between predicted and observed values.  The lines through each point give the standard deviation of the predictions for that point across all iterations.  There’s scatter here, but the model predicts low values for low observations, and high for high, so it’s a start…