Question

I have two files. One file has five observations of “Flux” across a three treatment experiment (treatments=A, B, C). In these three treatments temperatures have been manipulated. The observations of Flux are taken at five points in a 24 hr period. The second file (Temp) contains the temperatures for the three treatments across the 24 hour period.

I would like to use linear interpolation to predict what the Flux values will be at every hour during the 24 hour period. Note that the interpolation equations will be slightly different between the three treatments.

Can this be done in a loop so that the values of flux are estimated for each hour in the Temp.csv file? Then have the values integrated (summed) across the 24 hour period?

The files are available on dropbox here: Temp Data

This shows: the different slopes of the best fit linear relationships between flux and temperature across the three treatments:

#subset data in flux by treatment
fluxA<-flux[which(flux$Treatment=='A'),]
fluxB<-flux[which(flux$Treatment=='B'),]
fluxC<-flux[which(flux$Treatment=='C'),]

#Regression of Flux~Temperature
modelA<-lm (Flux~Temperature, data=fluxA)
summary (modelA)

modelB<-lm (Flux~Temperature, data=fluxB)
summary (modelB)

modelC<-lm (Flux~Temperature, data=fluxC)
summary (modelC)

#plot the regressions
plot (Flux~Temperature, data=fluxA,pch=16, xlim=c(0,28), ylim=c(0,20))
abline(modelA)

points(Flux~Temperature, data=fluxB,pch=16, col="orange")
abline(modelB, col="orange")

points(Flux~Temperature, data=fluxC,pch=16, col="red")
abline(modelC, col="red")
Was it helpful?

Solution

caldat <- read.csv(text="Treatment,Temperature,Flux
A,18.64,7.75
A,16.02,8.49
A,17.41,9.24
A,21.06,4.42
A,22.8,5.61
B,19.73,5.7
B,17.45,8.37
B,19.2,5.27
B,20.97,3.37
B,27.6,2.26
C,23.79,9.91
C,15.89,15.8
C,21.93,10.28
C,24.79,6.33
C,26.64,6.64
")

plot(Flux~Temperature, data=caldat, col=Treatment)
mod <- lm(Flux~Temperature*Treatment, data=caldat)
summary(mod)
points(rep(seq(16,28, length.out=1e3),3), 
       predict(mod, newdata=data.frame(Temperature=rep(seq(16,28, length.out=1e3),3),
                                       Treatment=rep(c("A", "B", "C"), each=1e3))),
       pch=".", col=rep(1:3, each=1e3))

enter image description here

You'll need to consider carefully if this is an appropriate and "good" model. Use standard regression diagnostics.

preddata <- read.csv(text="Time,A,B,C
100,17.8,21.64,23.04
200,17.5,21.3,22.7
300,17.23,21,22.39
400,16.92,20.67,22.08
500,16.47,20.3,21.74
600,15.78,19.75,21.24
700,15.19,19.14,20.63
800,14.58,18.47,20
900,14.22,17.99,19.49
1000,13.77,17.55,19.08
1100,13.39,17.02,18.62
1200,13.34,16.76,18.26
1300,13.17,16.62,18.05
1400,13.24,16.58,17.91
1500,13.31,16.63,17.86
1600,13.26,16.61,17.81
1700,13.12,16.57,17.75
1800,12.9,16.45,17.65
1900,12.74,16.32,17.54
2000,12.57,16.2,17.42
2100,12.36,16.04,17.28
2200,12.1,15.83,17.1
2300,11.79,15.57,16.88
2400,11.53,15.3,16.64
")

library(reshape2)
preddata <- melt(preddata, id="Time",
                 variable.name="Treatment", value.name="Temperature")
preddata$Flux <- predict(mod, newdata=preddata)

plot(Flux~Time, data=preddata, col=Treatment)

enter image description here

Sum the fluxes:

aggregate(Flux ~ Treatment, data=preddata, FUN=sum)
#  Treatment     Flux
#1         A 247.5572
#2         B 159.6803
#3         C 309.6186
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top