It's difficult to see where your error is when you provide no accompanying
data for us to test your code out on (X?  Y?  Stand?).

However, it looks like you are making your code more complex than it needs
to be.  Are you simply trying to a fit a separate linear regression to each
subset of your data identified by a unique standID?  If so, you can do so
like this:

# a fake version of your data
mydata <- data.frame(standID=rep(1:5, rep(20, 5)), X=rnorm(100),
Y=rnorm(100))

# split the data up by stand ID
standlist <- split(mydata, mydata$standID)

# fit a regression to each stand and save some results
standEst1 <- lapply(standlist, function(df) {
fit <- summary(lm(Y ~ X, data=df))
 co <- coef(fit)
intercept <- co[1, 1]
slope <- co[2, 1]
 pValue <- co[2, 4]
mse <- fit$sigma^2
rSquare <- fit$r.squared
 c(standID=df$standID[1], intercept=intercept, slope=slope, mse=mse,
rSquare=rSquare, pValue=pValue)
})

# convert the results to a data frame
standEstimates <- data.frame(do.call(rbind, standEst1))

Jean



On Thu, May 16, 2013 at 5:17 AM, rishard <or...@uclive.ac.nz> wrote:

> Hey I'm not really sure what I should put on here, but I am having trouble
> with my R code.  I am trying to get the p-values, R^2s etc for a number of
> different groups of variables that are all in one dataset.
>
> This is the code:
>
> #Stand counter
> st<-1
> #Collections
> stands<-numeric(67)
> slopes<-numeric(67)
> intercepts<-numeric(67)
> mses<-numeric(67)
> rsquares<-numeric(67)
> pValues<-numeric(67)
> #Start lists for X and Y values within each stand
> xi<-numeric(0)
> yi<-numeric(0)
> #Set the first element to the starting X and Y values
> xi[1]=X[1]
> yi[1]=Y[1]
> #Start looping working through your data, record by record
> for (i in 2:length(X)) {
>   #If you are in the same stand as on the last record, continue to
>   #collect X and Y values
>   if(Stand[i]==Stand[i-1]) {
>     xi=cbind(xi,X[i])
>     yi=cbind(yi,Y[i])
>   } else {
>     #If a new stand is encountered make your linear model and
>     #collect statistics
>     model<-lm(yi~xi)
>     stands[st]<-Stand[i-1]
>     intercepts[st]<-model$coefficients[1]
>     slopes[st]<-model$coefficients[2]
>     mses[st]<-sum(resid(model)^2)/(length(xi)-2)
>     ssr<-var(yi)*(length(xi)-1)-sum(resid(model)^2)
>     rsquares[st]<-ssr/(var(yi)*(length(xi)-1))
>     fRatio<-ssr/mses[st]
>     pValues[st]<-1-pf(fRatio,1,length(xi)-2)
>     #Increment the stand number, zero the within stand collections,
>     #and start again
>     st<-st+1
>     xi<-numeric(0)
>     yi<-numeric(0)
>     xi[1]=X[i]
>     yi[1]=Y[i]
>   }
> }
> #Make your data set
>
> standEstimates<-data.frame(standID=stands,intercept=intercepts,slop=slopes,mse=mses,rSquare=rsquares,pValue=pValues)
>
> The standEstimate outputs look like this:
>
> standID intercept       slop    mse     rSquare pValue
> 1       6833    319.2470        NA      0       NA      NA
> 2       756     708.7470        NA      0       NA      NA
> 3       795             508.2290        NA      0       NA      NA
> 4       1249    503.1460        NA      0       NA      NA
> 5       1331    703.0620        NA      0       NA      NA
> 6       1417    747.7620        NA      0       NA      NA
> 7       4715    549.3400        NA      0       NA      NA
> 8       4850    603.9940        NA      0       NA      NA
> 9       2105    573.3270        NA      0       NA      NA
> Etc etc
>
> and I get these warnings:
>
> 1: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 2: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>   number of items to replace is not a multiple of replacement length
> 3: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 4: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>   number of items to replace is not a multiple of replacement length
> 5: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 6: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>   number of items to replace is not a multiple of replacement length
> 7: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 8: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>   number of items to replace is not a multiple of replacement length
> 9: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 10: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>   number of items to replace is not a multiple of replacement length
> 11: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 12: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
>   number of items to replace is not a multiple of replacement length
> 13: In rsquares[st] <- ssr/(var(yi) * (length(xi) - 1)) :
>   number of items to replace is not a multiple of replacement length
> 14: In pValues[st] <- 1 - pf(fRatio, 1, length(xi) - 2) :
> etc etc
>
> Sorry if I haven't set this post out right, or haven't provided enough
> information.  But can anyone see why it is not giving me any returns for
> R^2
> etc?
>
>
>
> --
> View this message in context:
> http://r.789695.n4.nabble.com/R-looping-help-tp4667180.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@r-project.org 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.
>

        [[alternative HTML version deleted]]

______________________________________________
R-help@r-project.org 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.

Reply via email to