\[ \]

1. Data Step

There are 5 spots which have been selected.

\[ \]

1.1 SK-Old Aspen

Lat: 53.6289 / Long: -106.1978

Time range: 2003 – 2008

1.2 SK-Old Jack Pine

Lat: 53.9163 / Long: -104.6920

Time range: 2003 – 2008

1.3 SK-Southern Old Black Spruce

Lat: 53.9872 / Long: -105.1178

Time range: 2003 – 2008

1.4 SK-1975 (Young) Jack Pine

Lat: 53.8758 / Long: -104.6453

Time range: 2004 – 2006

1.5 QC-2000 Harvested Black Spruce/Jack Pine (HBS00)

Lat:49.26708 / Long:-74.03653

Time range: 2003 – 2010

1.6 Corresponding NDVI Data

\[ \]

2.Modelling Step

2.1 Separate Models for 5 spots

2.1.1 Old Aspen

NDVI <- OldAspen_NDVI
Indi0 <- c(rep("W",5),rep("S",15),rep("W",3)) 
Indi <- rep(Indi0,6)
Time <- rep(seq(1,366,by=16),6)

OldAspen.m <- gam(OldAspen_16day ~ s(Time, bs="cc", k=12) + NDVI + Indi)
gam.check(OldAspen.m)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradiant at convergence was 4.033404e-05 .
## The Hessian was positive definite.
## The estimated model rank was 13 (maximum possible: 13)
## Model rank =  13 / 13 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##            k'   edf k-index p-value
## s(Time) 10.00  9.39    1.20    0.99
#summary(OldAspen.m)
#OldAspen_var <- var(OldAspen.m$residuals)
OldAspen_var <- var(residuals(OldAspen.m))
#plot(OldAspen.m,residuals = TRUE,pch=19)
#plot(OldAspen.m$fitted.values,pch=19)
#lines(OldAspen.m$fitted.values)

2.1.2 Old Jack Pine

NDVI <- OldJackPine_NDVI
Indi0 <- c(rep("W",5),rep("S",15),rep("W",3)) 
Indi <- rep(Indi0,6)
Time <- rep(seq(1,366,by=16),6)

OldJackPine.m <- gam(OldJackPine_16day ~ s(Time, bs="cc", k=12) + NDVI + Indi)
gam.check(OldJackPine.m)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 6 iterations.
## The RMS GCV score gradiant at convergence was 3.086343e-07 .
## The Hessian was positive definite.
## The estimated model rank was 13 (maximum possible: 13)
## Model rank =  13 / 13 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##            k'   edf k-index p-value
## s(Time) 10.00  8.05    1.13    0.94
#summary(OldJackPine.m)
OldJackPine_var <- var(residuals(OldJackPine.m))

2.1.3 Southern Old Black Spruce

NDVI <- Southern_NDVI
Indi0 <- c(rep("W",5),rep("S",15),rep("W",3)) 
Indi <- rep(Indi0,6)
Time <- rep(seq(1,366,by=16),6)

Southern.m <- gam(Southern_16day ~ s(Time, bs="cc", k=12)+NDVI +Indi)
gam.check(Southern.m)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 6 iterations.
## The RMS GCV score gradiant at convergence was 4.179887e-07 .
## The Hessian was positive definite.
## The estimated model rank was 13 (maximum possible: 13)
## Model rank =  13 / 13 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##            k'   edf k-index p-value
## s(Time) 10.00  7.83    1.17    0.98
#summary(Southern.m)
Southern_var <- var(residuals(Southern.m))

2.1.4 Young Jack Pine

NDVI <- YoungJackPine_NDVI 
Indi0 <- c(rep("W",5),rep("S",15),rep("W",3)) 
Indi <- rep(Indi0,3)
Time <- rep(seq(1,366,by=16),3)

YoungJackPine.m <- gam(YoungJackPine_16day ~ s(Time, bs="cc", k=6) + NDVI + Indi)
gam.check(YoungJackPine.m)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradiant at convergence was 8.30027e-07 .
## The Hessian was positive definite.
## The estimated model rank was 7 (maximum possible: 7)
## Model rank =  7 / 7 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##           k'  edf k-index p-value
## s(Time) 4.00 3.59    1.07    0.67
#summary(YoungJackPine.m)
#YoungJackPine_var <- var(YoungJackPine.m$residuals)
YoungJackPine_var <- var(residuals(YoungJackPine.m))

2.1.5 Harvested Black Spruce/Jack Pine

NDVI <- QC2000_NDVI 
Indi0 <- c(rep("W",5),rep("S",15),rep("W",3)) 
Indi <- rep(Indi0,8)
Time <- rep(seq(1,366,by=16),8)

QC2000.m <- gam(QC2000_16day ~ s(Time, bs="cc", k=12) + NDVI + Indi)
gam.check(QC2000.m)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 10 iterations.
## The RMS GCV score gradiant at convergence was 1.120602e-07 .
## The Hessian was positive definite.
## The estimated model rank was 13 (maximum possible: 13)
## Model rank =  13 / 13 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'      edf  k-index p-value
## s(Time) 1.00e+01 9.19e-09 1.16e+00    0.99
#summary(QC2000.m)
#QC2000_var <- var(QC2000.m$residuals)
QC2000_var <- var(residuals(QC2000.m))

2.2 Combine 5 spots

FlxTwr <- c(OldAspen_16day,OldJackPine_16day,Southern_16day,YoungJackPine_16day,QC2000_16day)
ndvi <- c(OldAspen_NDVI,OldJackPine_NDVI,Southern_NDVI,YoungJackPine_NDVI,QC2000_NDVI)

time <- rep(seq(1,365,by=16),29)
Indi0 <- c(rep("W",5),rep("S",15),rep("W",3)) 
indi <- rep(Indi0,29)

spe <- c(rep("OldAspen",length(OldAspen_NDVI)),rep("OldJackPine",length(OldJackPine_NDVI)),
           rep("BlackSpr",length(Southern_NDVI)),rep("YoungJackPine",length(YoungJackPine_NDVI)),rep("BlackSpr",length(QC2000_NDVI)) )

Weight <- c(rep(1/OldAspen_var,138),rep(1/OldJackPine_var,138),rep(1/Southern_var,138),rep(1/YoungJackPine_var,69),rep(1/QC2000_var,184))
w <- Weight/mean(Weight)

D <- data.frame(Flx=FlxTwr,Time=time,NDVI=ndvi,Spe=spe,W=w,Indi=indi)

rownames(D) <- paste0("row",1:nrow(D)) 
Mix <- gam(Flx ~ s(Time,bs="cc",k=23)+s(Time,bs="cc",k=23,by=Spe)+ NDVI + Indi ,weights = W,data = D)
Diff <- D[setdiff(rownames(D),names(predict(Mix))),]
Diff.row <- as.numeric(substr(rownames(Diff),4,6))

Mplot <- function(data,name){
  label <- c("Jan 1","Jan 17","Feb 2","Feb 18","Mar 6","Mar 22","April 7","April 23","May 9","May 25","Jun 10","Jun 26","July 12","July 28","Aug 13","Aug 29","Sep 14","Sep 30","Oct 16","Nov 1","Nov 17","Dec 3","Dec 19")
  monthplot(data,phase = seq_along(data)%%23+1,labels=label,ylab="FlxTwr")
  abline(v=9,lty="dashed")
  title(name)
}

Smooth and Registration

  daytime <- seq(1,365,by=16)
  OA <- matrix(OldAspen_16day,nrow=23)
  OAbasis <- create.bspline.basis(c(1,365),nbasis=13,norder=4,breaks=c(1,120,130,140,170,200,220,250,260,270,365))
  OAbasismat <- eval.basis(daytime, OAbasis)
  OAcoef <- solve(crossprod(OAbasismat,OAbasismat))%*%crossprod(OAbasismat,OA)
  everyday <- 1:365
  yhat <-  eval.basis(everyday,OAbasis)%*%OAcoef

  plot(yhat[,1],type="l",ylim=c(-6,6),main="Curve after smooth only",
  xlab="Day of the year")
for(i in 2:6){
  lines(yhat[,i],col=i)
}

Plot_1st_deri <- function(){
    OAbasismat_1st <- eval.basis(everyday,OAbasis,1)
    y_1st <- OAbasismat_1st%*%OAcoef
    for(i in 1:6){
      plot(yhat[,i],type="l",ylim=c(-6,3),main=paste("First derivative of Year",i))
      points(y_1st[,i]*10)
      abline(h=0,col=2)
    }
}

App <- function(TStars,Gpoints){
  approx(c(1,TStars,365),c(1,Gpoints,365),everyday)$y
}

Regi <- function(Ti,yi,basis=OAbasis){
  OAcoef <- solve(crossprod(OAbasismat,OAbasismat))%*%crossprod(OAbasismat,yi)
  eval.basis(Ti,basis)%*%OAcoef
}


#Gmat <- Gmat_3landmarks
Gmat0 <- read.csv("/Users/fhe7/Documents/School/Oct_15 Gammon/G3landmarks.csv",header=T)
Gmat <- Gmat0[,-1]
  TStars <- apply(Gmat,2,mean)
  Tmat <- matrix(NA,nrow=6,ncol=365)
  for(i in 1:6){
    Tmat[i,] <- App(TStars,Gmat[i,])
  }
  
  Regimat <- matrix(NA,nrow=6,ncol=365)
  for(i in 1:6){
    Regimat[i,] <- Regi(Tmat[i,],OA[,i])
  }
  
  plot(1:365,Regimat[1,], type="l", ylim=c(-6,6) ,main="Curve after registration at 3 landmarks",xlab="Day of the year")
  
  for(i in 2:6){
    lines(1:365,Regimat[i,])
  }
  abline(v=TStars,lty="dashed")

P1: Find curve for future year.

[1] Pick some of the T* (let’s say \(n\) of them, \(n = 10\) or more).

Evaluate each year at T*\(_1\) till T*\(_n\).

[2] Choose \(m\) landmarks.

[3] Each curve now has \(m+n\) values.

[4] Find the joint distribution.

[4-1] obs number \(= (n+m) x 6\), as for the unobserved year 7, sample from multivariate normal

[4-2] Raw points \(\rightarrow\) Regi curve (shift) \(\rightarrow\) Fit again (fewer knots?) \(\rightarrow\) using the coefs?

P2: Ambiguity/Bad landmarks.

P3: Distin MC VS. non-homogenious-Poisson Process