I have a problem which is extra irritating, as my colleague gets all done properly in an old standalone version of JMP, but I can't get it right in R :-(
I have a dataset with measurements per site per year:
SITE Year N
1 8 1990 7.33
2 8 1991 7.30
3 8 1992 7.00
4 8 1993 7.49
5 8 1994 7.91
6 8 1995 7.76
etc.
The file can be downloaded from https://dl.dropboxusercontent.com/u/14125905/N_test.csv
# Import data
> Nind_SiteValues <-
> read.table("N_test.csv",
> header=TRUE, sep=";", na.strings="NA", dec=".", strip.white=TRUE)
I want to see if there is a trend in R and make a graph of N per year. As sites differ per year, I have to add them as random factor:
# Use lme with SITE as random factor
> library(nlme)
> N_trend <- lme(N ~ Year, random=~1 | factor(SITE), data=Nind_SiteValues)
> summary(N_trend)
Linear mixed-effects model fit by REML
Data: Nind_SiteValues
AIC BIC logLik
19163.11 19191.49 -9577.553
Random effects:
Formula: ~1 | factor(SITE)
(Intercept) Residual
StdDev: 1.381553 0.5550342
Fixed effects: N ~ Year
Value Std.Error DF t-value p-value
(Intercept) -18.266466 2.5542375 7643 -7.151436 0
Year 0.012033 0.0012736 7643 9.447777 0
Correlation:
(Intr)
Year -1
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-5.712515639 -0.477779874 -0.005360463 0.430617650 10.276334141
Number of Observations: 8922
Number of Groups: 1278
So this gives me a slope of 0.012033 which is exactly the same as my colleage in JMP. So far so good. When I try to get the values per year for the graph I thought I had to use:
# Calculate values per year for graph
> N_year <- lme(N ~ factor(Year), random=~1 | factor(SITE), data=Nind_SiteValues)
> summary(N_year)
Linear mixed-effects model fit by REML
Data: Nind_SiteValues
AIC BIC logLik
19021.52 19220.13 -9482.76
Random effects:
Formula: ~1 | factor(SITE)
(Intercept) Residual
StdDev: 1.383604 0.5450457
Fixed effects: N ~ factor(Year)
Value Std.Error DF t-value p-value
(Intercept) 5.696596 0.07563604 7619 75.31589 0.0000
factor(Year)1991 -0.084066 0.07765718 7619 -1.08253 0.2791
factor(Year)1992 0.068109 0.07379237 7619 0.92298 0.3560
factor(Year)1993 -0.036885 0.07388291 7619 -0.49924 0.6176
factor(Year)1994 -0.001379 0.07392474 7619 -0.01865 0.9851
factor(Year)1995 0.085483 0.07394456 7619 1.15604 0.2477
factor(Year)1996 -0.158037 0.07356505 7619 -2.14826 0.0317
factor(Year)1997 0.192585 0.07331735 7619 2.62674 0.0086
factor(Year)1998 0.140250 0.07293321 7619 1.92299 0.0545
factor(Year)1999 0.096754 0.07324577 7619 1.32094 0.1866
factor(Year)2000 0.149507 0.07327045 7619 2.04048 0.0413
factor(Year)2001 0.114410 0.07331167 7619 1.56060 0.1187
factor(Year)2002 0.236566 0.07328951 7619 3.22783 0.0013
factor(Year)2003 0.147536 0.07157412 7619 2.06130 0.0393
factor(Year)2004 -0.063435 0.07145848 7619 -0.88772 0.3747
factor(Year)2005 0.167900 0.07142229 7619 2.35080 0.0188
factor(Year)2006 0.315268 0.07141166 7619 4.41480 0.0000
factor(Year)2007 0.332034 0.07155791 7619 4.64008 0.0000
factor(Year)2008 0.098167 0.07165542 7619 1.36999 0.1707
factor(Year)2009 0.161026 0.07168395 7619 2.24633 0.0247
factor(Year)2010 0.199272 0.07155167 7619 2.78501 0.0054
factor(Year)2011 0.428803 0.07149717 7619 5.99748 0.0000
factor(Year)2012 0.356957 0.07157344 7619 4.98728 0.0000
factor(Year)2013 0.263030 0.07153716 7619 3.67682 0.0002
factor(Year)2014 0.163145 0.07184643 7619 2.27074 0.0232
factor(Year)2015 0.220120 0.07214547 7619 3.05105 0.0023
The graph of slope and year values looks completely wrong (upper graph in https://dl.dropboxusercontent.com/u/14125905/Graphs.png)
The graph of the JMP version of my colleague looks like the lower graph (so my slope, his year values).
Obviously there is something wrong in the way I calculate the year values in R, but I can't figure out what. Also strange that there is a perfect match for the slope, not for the values. Any help or guidance would be appreciated.
