> load(file.choose()) #grab new Tf3 object from course webpage

> Tf3$Date = as.Date(Tf3$DATE,format="%m/%d/%Y")  #making day a 'Date' object

> Tf4 = na.exclude(Tf3)        #remove missing values

> library(nlme)

 

 

Adding a spatial component to the data, as is, was possible because the dataset had the UTM coordinates being repeated over the various dates. So, one transect was used instead.

 

1.       Adding a spatial component to temporal dataset (using Date as the random effect variable):

>mod8 = lme(fixed=minT~minSYN,random=~1|DATE,data=Tf4,correlation=corSpher(c(1,0.9),form=~UTME_CORR+UTMN_CORR|DATE,nugget=T))

 

Warning message:

In Initialize.corSpher(X [[2L]] ...):

  Initial value for range less than minimum distance. Setting it to 1.1 * min (distance)

 

> summary (mod8)

 

Linear mixed-effects model fit by REML

 

Data: Tf4

       AIC                  BIC                         logLik

  7249.719            7286.292              -3618.859c

 

Random effects:

 Formula: ~1 | DATE

                                (Intercept)          Residual

StdDev:                 1.798772             1.141832

 

Correlation Structure:   Spherical spatial correlation

 Formula: ~UTME_CORR + UTMN_CORR | DATE

 Parameter estimate(s):

      range              nugget

714.7034439       0.0997957

Fixed effects: minT ~ minSYN

                                Value                    Std.Error              DF           t-value                                 p-value

(Intercept)                          1.719643              0.11771657          2982       14.60834               0

minSYN                                 0.878300             0.01290723          2982       68.04713              0

 

Correlation:

       (Intr)

minSYN -0.202

 

Standardized Within-Group Residuals:

                 Min                       Q1                          Med                       Q3                          Max

-3.9235741          -0.2254490          0.1928420            0.5448933            3.1102055

 

Number of Observations: 3282

Number of Groups: 299

 

2.       Comparing the model with a correlation structure to one without a correlation structure:

 

> mod7 = update (mod8, correlation=NULL)

 

> summary (mod7)

Linear mixed-effects model fit by REML

 Data: Tf4

       AIC      BIC    logLik

  8476.279 8500.661 -4234.139

 

Random effects:

 Formula: ~1 | DATE

                                 (Intercept)         Residual

StdDev:                2.042481              0.715257

 

Fixed effects: minT ~ minSYN

                                Value                    Std.Error              DF           t-value                                 p-value

(Intercept)          1.8333195            0.1211843            2982       15.12836              0

minSYN                0.9050649            0.0128835            2982       70.24995              0

 

Correlation:

       (Intr)

minSYN -0.198

 

Standardized Within-Group Residuals:

                 Min                       Q1                          Med                       Q3                          Max

-5.27767740        -0.47945199        0.05280634          0.50181355          3.67454470

 

Number of Observations: 3282

Number of Groups: 299

 

> AIC(mod7,mod8)

                                df            AIC

mod7                    4              8476.279

mod8                    6              7249.719  

 

> anova(mod7,mod8)

 

Model   df            AIC                         BIC                         logLik                    Test  L.Ratio       p-value

mod7    1             4              8476.279              8500.661              -4234.139                       

mod8     2             6              7249.719 7           286.292                 -3618.859             1 vs 2 1230.560 <.0001

 

 

 

3.       Comparing Variograms:

 

> plot (Variogram(mod7,form=~UTME_CORR+UTMN_CORR))

 

Mod7: Variogram without a correlation structure

 

> plot (Variogram(mod8,form=~UTME_CORR+UTMN_CORR))

 

 

 

 

 

 

 

 

 

 

 

 

4.       Adding a Exponential correlation to see if this would improve the model:

 

> mod9 = update (mod8, correlation=corExp(c (1, 0.9), form=~UTME_CORR+UTMN_CORR|DATE, nugget=T))

 

> summary(mod9)

Linear mixed-effects model fit by REML

 Data: Tf4

                 AIC                        BIC                         logLik

                8480.279              8516.852              -4234.139

 

Random effects:

 Formula: ~1 | DATE

        (Intercept) Residual

StdDev:    2.042481 0.715257

 

Correlation Structure: Exponential spatial correlation

 Formula: ~UTME_CORR + UTMN_CORR | DATE

 

 Parameter estimate(s):

 Range nugget

   1.0    0.9

 

Fixed effects: minT ~ minSYN

                                Value                    Std.Error              DF           t-value                                 p-value

(Intercept)          1.8333195            0.1211843            2982       15.12836              0

minSYN                0.9050649            0.0128835            2982       70.24995              0

 

 Correlation:        (Intr) minSYN -0.198

 

Standardized Within-Group Residuals:

                Min                        Q1                          Med                       Q3                          Max

-5.27767740        -0.47945199        0.05280634          0.50181355          3.67454470

 

Number of Observations: 3282

Number of Groups: 299

 

> AIC(mod8,mod9)

                 df           AIC

mod8     6             7249.719

mod9    6              8480.279

 

 

> anova (mod8,mod9)

                Model   df            AIC                         BIC                         logLik

mod8     1              6             7249.719              7286.292              -3618.859

mod9     2              6             8480.279             8516.852              -4234.139