RDX2 X  errorsINxGP source" +function(n=1000, xshift=1.25, a=2.5, b=1.5, + gpSDx=2, yerr=0.5, sdRatio=1.5){ I ## xshift is the difference in mean value of x, between ctl and TRT # xraw <- rnorm(2*n, 12.5, gpSDx) J df <- data.frame(x=rnorm(2*n, 12.5, gpSDx) + rep(c(0,xshift), c(n,n)), 8 trt = rep(c("CTL", "TRT"), c(n,n))) , df$y <- a + b*df$x + rnorm(2*n, sd=yerr) C ## b is the slope of the regression line of y on the "true" x < ## Now add random error to each of the $x$-values, thus: E ## gpSDx is the conditional SD of x, given the treatment groups 6 df$xWITHerr <- df$x + rnorm(2*n, sd=sdRatio*gpSDx) @ ## sdRatio is the ratio of the measurement error SD to SDx $ b <- lm(y ~ trt+x, data=df)$coef + B <- lm(y ~ trt+xWITHerr, data=df)$coef 7 atten <- 1/(1+sdRatio^2) # Theoretical attenuation ? cat("\n", "Attenuation, simulated model =", round(atten,3)) I cat("\n", "Slope Estimates - On x:", round(b[3],3), " On xWITHerr:",  round(B[3],3), "\n") A cat("\n", "Estimated group effect - Reg on x", round(b[2],3), 1 " Reg on xWITHerr", round(B[2],3), "\n") 6 plt <- xyplot(y ~ x+xWITHerr, data=df, groups=trt, * panel=function(x,y,...){ 5 panel.superpose(x,y,type="p",...) 0 subs <- list(...)$subscripts 1 fac <- list(...)$groups[subs] ) gps.lm <- lm(y~fac+x) ) hat <- fitted(gps.lm) 9 panel.superpose(x,y=hat,type="r",...)  }, * auto.key=list(columns=2)  )  print(plt)  invisible(df)  }þý n@@ xshift?ô a@ b?ø gpSDx@ yerr?à$ sdRatio?øþ srcref  # srcfileò encoding native.enc timestampAÑwç€ class POSIXt POSIXctþ filename ~/r-book/ed2/xtras/errorsINxGP.R wd /Users/johnm/r/lundþþÿ srcfileþÿ srcrefþ 8 ÿ ÿÿ srcrefþ , ÿ ÿÿ srcrefþ   6 ÿ ÿÿ srcrefþ $ ÿ ÿÿ srcrefþ + ÿ ÿÿ srcrefþ  ÿ ÿÿ srcrefþ ? ÿ ÿÿ srcrefþ  ÿ ÿÿ srcrefþ 1 ÿ ÿÿ srcrefþ   ÿ ÿÿ srcrefþ !! ÿ ÿÿ srcrefþ "" ÿ ÿÿ srcrefþþ { <- xraw rnorm *@ÿþ@)ÿþþÿ df data.frame x +ÿÿ@ÿþ@)ÿþ rep cÿþÿÿÿþþþ trtÿÿ CTL TRTþÿÿÿþþþþÿ $ÿ yþÿÿÿÿÿÿÿÿþþþÿÿ@ÿþ sdÿþþþÿÿÿ xWITHerrþÿÿÿÿþÿÿ@ÿþ ÿÿ ÿÿþþþþÿÿÿ lm ~ÿÿÿÿþþ dataÿþ coefþþÿ Bÿ"ÿ#ÿÿÿÿ!ÿþþ$ÿÿþ%ÿþþÿ atten /?ð (ÿ?ð ^ ÿ@þþþþþ cat   Attenuation, simulated model = round'ÿ@þþ+ÿ   Slope Estimates - On x:,ÿ [ÿ@þ@þ  On xWITHerr:,ÿ-ÿ&ÿ@þ@þ  þ+ÿ   !Estimated group effect - Reg on x,ÿ-ÿÿ@þ@þ  Reg on xWITHerr,ÿ-ÿ&ÿ@þ@þ  þÿ plt xyplot#ÿÿÿÿ!ÿþþ$ÿÿ groupsÿ panel function$ÿû$ÿû ...ûþ ÿ 5 ÿ ÿÿ srcrefþ 0 ÿ ÿÿ srcrefþ 1 ÿ ÿÿ srcrefþ ) ÿ ÿÿ srcrefþ ) ÿ ÿÿ srcrefþ 9 ÿ ÿÿ srcrefþþÿ panel.superposeÿÿ type p3ÿþÿ subsÿ list3ÿþ subscriptsþþÿ fac-ÿÿ7ÿ3ÿþ0ÿþ6ÿþþÿ gps.lm"ÿ#ÿÿÿ9ÿÿþþþþÿ hat fitted:ÿþþ4ÿÿÿ;ÿ5ÿ r3ÿþþ function(x,y,...){ 5 panel.superpose(x,y,type="p",...) 0 subs <- list(...)$subscripts 1 fac <- list(...)$groups[subs] ) gps.lm <- lm(y~fac+x) ) hat <- fitted(gps.lm) 9 panel.superpose(x,y=hat,type="r",...)  }þ auto.key7ÿ columns@þþþ print.ÿþ invisibleÿþþþ