In R, the output of my linear model shows a positive correlation but my ggplot graph indicates a negative correlation?

I suppose this is because you look at the raw data you fed into the model, not the model predictions. In the plot, you don't "isolate" a single predictor. You look at the result of all predictors doing something to the response variable. I suppose the effect of this predictor is "overshadowed" by the effects of the other predictors.

To take a look at the effect of solely one predictor, you need to predict new values from the model while fixing all other predictors. You can try something along the lines of:

  preds <- predict(yourmodel, newdata = data.frame(
    "Sycamore_biomass" = 0:25,
    "Distance_to_road" = mean(egglay_date$Distance_to_road),
    "Distance_to_light" = mean(egglay_date$Distance_to_light),
    "Anthropogenic_cover" = mean(egglay_date$Anthropogenic_cover),
    "Canopy_cover" = mean(egglay_date$Canopy_cover)))
  
  new_data <- data.frame(
    "Sycamore_biomass" = 0:25,
    "First_egg" = preds)
  
  ggplot(new_data, aes(x=Sycamore_biomass, y=First_egg)) +
    geom_point(shape=19, alpha=1/4) +
    geom_smooth(method=lm)

This should give you the predictions of your model when only considering the effect of the one predictor.

The answer to your question is quite simple (but I understand why it may seems complex at first).

First off, your model indicates a positive relationship because you have included all your other variables. Keep in mind, your best fit line through your data here is when you take all your data points, and fit a line to make the sum of residuals = 0. Note: this is not the same as sum of residuals squared.

Since you didn't provide your data (please do on future posts, or at least, something to work with), I will illustrate my point with the data(mtcars) built into R

data("mtcars")

df <- mtcars

This dataset has many variables, to see them all, just type names(df)

Lets just work with three of them to see if miles per gallon (mpg) is explained by:

   1) cyl : # of cylinders
   2) hp : horse power
   3) drat : rear axle ratio

attach(df)
model <-lm(mpg~cyl+hp+drat)
summary(model)

Let's say, I just want to plot the relationship between cylinders and mpg (for you, it would be sycamour biomass and bird lay). Here, from our model summary, we see that our relationship is negative (negative estimate, aka, coefficient), and that the intercept is at 22.5.

So I do what you just did and just plot mpg~cly (without considering my other variables)

plot(mpg~cyl, pch=15, col="blue",cex=2, cex.axis=2, ylab="MPG", xlab="Number of Cylinders", cex.lab=1.5)
abline(lm(mpg~cyl),lwd=2,col="red")

First off, we see that the y intercept is not 22.5, but rather above 25.

If I were to do the math from first model, if I had 4 cylinders, I should predict:

22.51406 + (4 * -1.3606) = 17.07

So lets see if our prediction is correct on our graph

Definitely not.

So lets run a new model (which you need to do), where we model just mpg~cly

reduced_model <- lm(mpg~cyl)

summary(reduced_model)

See how the intercept and coefficent (estimates) changed? Yours will too when you run a reduced model. Lets see if the plots now make sense following the same steps as above with predicting 4 cylinders

37.8846  + (4 * -2.8758 ) # 26.38
plot(mpg~cyl, pch=15, col="blue",cex=2, cex.axis=2, ylab="MPG", xlab="Number of Cylinders", cex.lab=1.5)
abline(lm(mpg~cyl),lwd=2,col="red")
abline(h=26.38,v=4,lwd=2, col="green")

Looks like everything checks out.

Summary: You need to run a simple model with just your two variables of interest if you want to correctly understand your plot