[R] optimize calculations

From: Bart Joosen <bartjoosen_at_hotmail.com>
Date: Sun 01 Apr 2007 - 19:02:42 GMT


Hi,

I have a model, which has a logaritmic response, but there is an offset for this response. I use optimize to find a value to minimise the residuals of the model:

f2 <- function(x,df) sum(summary(lm(log((y-x)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df))$residuals^2) start <- optimize(f2,c(0,min(y)),df=df)$minimum mod <- lm(log((y-start)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df) As I have already some data, I have build my model, but I would like to simulate the effect of adding more data points on the confidence interval for the predicted values.

So my code was:

forcast <- function(mod2,start, df,stud.temp = 25,stud.end=52,sims=100,Temp=c(25,40), Mnd=c(0.5,1,1.5),sd=0.005){

    df <- as.data.frame(expand.grid(Mnd,Temp))     names(df) <- c("Mnd","Temp")
    df$Pred <- exp(predict(mod2,list(Temp= df$Temp, Mnd=df$Mnd)))*df$Mnd+start     res.int <- c(1:sims)
    f2 <- function(x,df) sum(summary(lm(log((y-x)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df))$residuals^2)     for (i in 1:sims){

        y <- rnorm(nrow(df),df$Pred,sd)
        start <- optimize(f2,c(0,min(y)),df=df)$minimum
        mod <- lm(log((y-start)/Mnd) ~ I(1e+05/(8.617*(Temp + 273.16))), df)
        pr <- (exp(predict(mod,list(Temp= stud.temp, Mnd=stud.end),interval="confidence"))*stud.end+start)
        res.int[i] <- (pr[,3]-pr[,2])/2

    }
    cat(Mnd,round(mean(res.int),3),"\n") }

Where one can input the base model, which is used for the prediction of the values for the generated expand.grid dataframe. For this predicted values, we use the rnorm function to simulate the error, and then the model is fitted for this values, and the width of the confidence interval at the given temperature and time point is calculated.

As I would like to investigate the effect of the number of months used on the width of the confidence interval, I would like to speed up the coding, so I can use more simulations to investigate the effect.

Thanks in advance

Bart

        [[alternative HTML version deleted]]



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon Apr 02 05:09:21 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 01 Apr 2007 - 19:30:37 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.