2
$\begingroup$

I'm trying to fit a logistic curve to cumulative data, derived from satellite imagery. Previously, I have point observation data which were either 0s or 1s. Os being 'forest' and 1s being 'non-forest'. These point observations existed for multiple images/dates. So I had one csv file with 'observation date' in once column and 'state' in another. Weights were also included, but these aren't relevant here.

This was the code I used:

# Read in data
data=read.csv("C:\\Users\\Owner\\Desktop\\BV\\Deforestation_Analysis\\Ambaro_Ambanja\\Analysis\\CDM_Table_Final.csv")

head(data)
attach(data)
Time<-strptime(Observation_Date, "%m/%d/%Y")
Time2<-as.Date(Time)
Time2
# Set project start date to zero
Zero<-as.Date("2015/7/1")
temp.time<-as.numeric(Zero-Time2)
T.Time<-temp.time*(-1)
T.Time
data2<-cbind(data,T.Time)

# Fit the model and summarise
model.glm=glm(formula=State~T.Time,data=data2,weights=Weight,family=binomial(link=logit))
summary(model.glm)

Now, rather than having multiple 0s and 1s for each date, I have simply number of non-forest pixels (and no weight). Here is the data I have:

enter image description here

The 'State' field is literally the portion of non-forest pixels (18.6% of pixels were non-forest in 2014)

I tried runnning the following adaption of the original script:

# Read in data
data=read.csv("C:\\Users\\Leah\\Desktop\\BV\\AAB\\Geospatial\\Deforestation_Analysis\\Analysis\\For_R.csv")

head(data)
attach(data)
Time<-strptime(Observation_Date, "%d/%m/%Y")
Time2<-as.Date(Time)
Time2
# Set project start date to zero
Zero<-as.Date("2015/7/1")
temp.time<-as.numeric(Zero-Time2)
T.Time<-temp.time*(-1)
T.Time
data2<-cbind(data,T.Time)

# Fit the model and summarise
model.glm=glm(formula=State~T.Time,data=data2,family=binomial(link=logit))
summary(model.glm)

I fully expected it to fail, because the data is no longer binomial. But while it did throw up a warning ('In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!'), the coefficients it spat out formed a curve which looked perfect. But I wasn't overly confident in this hash.

I've read that you can still use the binomial glm family if you feed R with a table containing successes (non-forest) and failures. So I came up with this adapted data:

enter image description here

and the following adapted script:

# Read in data
data=read.csv("C:\\Users\\Leah\\Desktop\\BV\\AAB\\Geospatial\\Deforestation_Analysis\\Analysis\\For_R_Final.csv")

head(data)
attach(data)
Time<-strptime(Observation_Date, "%d/%m/%Y")
Time2<-as.Date(Time)
Time2
# Set project start date to zero
Zero<-as.Date("2015/7/1")
temp.time<-as.numeric(Zero-Time2)
T.Time<-temp.time*(-1)
T.Time
data2<-cbind(data,T.Time)

# Fit the model and summarise
model.glm=glm(formula=cbind(Deforested, Total-Deforested)~T.Time,data=data2,family=binomial(link=logit))
summary(model.glm)

It ran fine with no errors, but the trend it generated doesn't fit the data anywhere near as well as the hashed version:

Blue line is hased version; Grey line is adapted script; Red points are the data I'm using to fit the model

The blue line is hased version; grey line is adapted script; red points are the data I'm using to fit the model.

Why does the adapted version fit the point worse than the hashed version? Is R so clever that it just uses my fractional values how I want it to in the glm(family=binomial)?

Any advice greatly appreciated! Am not happy with how I got the blue trend and this work is very important to our study.

THANK YOU!!

$\endgroup$
  • 3
    $\begingroup$ A warning about something people often miss about cumulative data -- Cumulative data tends to be highly serially correlated. If the increments are independent, then cumulating many of them can yield correlations quite close to 1. $\endgroup$ – Glen_b Oct 7 '14 at 10:10

Your Answer

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