tag:blogger.com,1999:blog-73161587356028134862015-09-16T17:10:04.221-07:00Eye HealthEye Statshttp://www.blogger.com/profile/17844914115422103029noreply@blogger.comBlogger5125tag:blogger.com,1999:blog-7316158735602813486.post-43584543128586352902014-09-02T21:37:00.001-07:002014-09-02T21:37:20.375-07:00Plaquenil ScreeningPlaquenil (hydroxychloroquine) and its cousin chloroquine are useful drugs - particularly for arthritis and some skin conditions. Both drugs can cause significant damage to the retina which is usually irreversible. Doctors have known for decades about the importance of monitoring of all patients on long term plaquenil.<br /><br />The technology for screening for plaquenil toxicity has evolved over recent years. A nice review article can be found <a href="http://www.aao.org/publications/eyenet/201105/upload/CUComp-May-2011.pdf" target="_blank">here</a>. There has been much enthusiasm recently for high resolution OCT scanning for detection of early plaquenil toxicity and this has in many doctors' practices been at the expense of visual field testing.<br /><br />A case was recently seen where substantial and permanent visual loss resulted from plaquenil despite the regular OCT scans.<br /><br />Our approach is now, where possible, do some baseline tests before the first dose of plaquenil. This includes colour vision, central field and OCT. We have found Fundus Autofluorescence (FAF) to be more effective than OCT in the detection of damage so for patients with more than 3 years of plaquenil treatment we recommend annual visual fields, FAF and OCT.<br /><br />In cases where patients have been on Plaquenil and it is working for their joint symptoms but there is some suspicion of toxicity, we recommend the addition of a <a href="http://www.essendonretina.com.au/index.php/electrophysiology/more-information" target="_blank">multifocal ERG test</a> . This is a test of the function of the retina rather than its structure which is tested by OCT and FAF<br /><br />At our <a href="http://www.essendoneye.com.au/" target="_blank">Essendon rooms</a> we can offer all of these types of tests as well as testing for a large range of other eye diseasesEye Statshttp://www.blogger.com/profile/17844914115422103029noreply@blogger.com0tag:blogger.com,1999:blog-7316158735602813486.post-6420686273316194222013-02-20T00:19:00.000-08:002013-02-20T01:16:01.099-08:00Bolt vs Emerson and the pause - who is right. Using piece-wise regression to find an answerThere has been some debate about the interpretation of global temperature readings - specifically whether there has been a "pause" in the rate of "global warming" over the last 16 years. An recent episode of this is presented at <a href="http://joannenova.com.au/2013/02/the-emerson-v-bolt-argument-on-air-does-emerson-not-know-statistics/">JoNova</a>. The actual exchange is available following the links at the above page. Leaving aside the important point of whether a global temperature is a valid measure, the issue comes down to how to statistically pose the question. Fitting a regression line to a selected time period and inferring that there is no significant trend is not really valid. Of course, fitting a linear regression line to the entire temperature record can be done <br /><br /><div style="text-align: center;"><a href="http://2.bp.blogspot.com/-lYk2G3binCk/USSLGGFu12I/AAAAAAAAAB4/satN8UsQiUk/s1600/Rplot02.png" imageanchor="1"><img border="0" src="http://2.bp.blogspot.com/-lYk2G3binCk/USSLGGFu12I/AAAAAAAAAB4/satN8UsQiUk/s320/Rplot02.png" /></a> </div>and you get a reasonable regression line. This model however assumes a stationary process over the period - no hockey stick and assumes a linearly increasing "human influence" over the entire period. It is actually a fairly good fit to the data.<br /><br /> The regression data from R are<br /><table border="1"><tbody><tr><td colspan="2">Coefficients</td></tr><tr><td>Intercept</td><td>Gradient</td></tr><tr><td>5.769e-02</td><td>1.272e-05</td></tr></tbody></table><br />The correlation coefficient is 0.744<br /><br />Of course, it looks much more dramatic if you restrict the data to post 1974. A fit to that data of a linear model gives<br /><table border="1"><tbody><tr><td colspan="2">Coefficients</td></tr><tr><td>Intercept</td><td>Gradient</td></tr><tr><td>-1.638e-01</td><td>4.641e-05</td></tr></tbody></table><br /> The correlation coefficient is 0.821<br /><br /> It correlates better and the gradient is 3.8 times as steep. It also corresponds to a plausible period of increased "human influence". A problem is that if you allow a change in the regression line at that time point, you must allow other changes as well and this poses the question of whether there has been another change and if so to what ?<br /><br />It is true that "eyeballing" the data there seems to be a flattening off in later years<br /><br />This graph looks at the data post 1974 and compares two models : one with a breakpoint and one without<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-yq64D6SvEvU/USSR0jThBQI/AAAAAAAAACI/-Hj50xesYcs/s1600/Rplot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="286" src="http://3.bp.blogspot.com/-yq64D6SvEvU/USSR0jThBQI/AAAAAAAAACI/-Hj50xesYcs/s320/Rplot.png" width="320" /></a></div><div class="separator" style="clear: both; text-align: left;">Now, comparing model fits like this is not easy. The likelihood function of the data fit needs to be adjusted for the number of parameters. The single line model has only 2 parameters - slope and intercept but the 2 line model has 5 parameters - the slope and intercept of both lines and the break point position. One way to compare the models is to use the AIC (Akaike Information Criterion) - the lower the better. In this case, the AIC for the single line was -568 and the two line model was -602. Not much difference but some evidence for a breakpoint in September 2005</div><div class="separator" style="clear: both; text-align: left;"><br /></div><br />Eye Statshttp://www.blogger.com/profile/17844914115422103029noreply@blogger.com0tag:blogger.com,1999:blog-7316158735602813486.post-6011349471153130872012-09-01T17:59:00.001-07:002012-09-01T17:59:15.836-07:00AstigmatismAstigmatism is a very difficult concept to explain to students and patients. There are no good webpage references on the combination of different spherocylinder lenses. Understanding how they add is however important in understanding how the Jackson Cross Cylinder works and also in understanding how astigmatic changes come about with surgery.<br />At <a href="http://nweye.com.au/" target="_blank">North West Eye Specialists</a> we have the <a href="http://www.hicsoap.com/">Holladay software</a> for analysing surgical astigmatism which is allowing each of our doctors to form their own database parameters for the optimized management of astigmatism.<br />As a teaching tool, i have developed an <a href="http://bees.net.au/eyehealth/spherocyl.htm">online astigmatism calculator</a> – it can add and subtract spherocylindrical lenses. It is still in beta form and not 100% accurate but will help in the understanding of astigmatism. Please feel free to leave suggestions or drop me a line for the source codeEye Statshttp://www.blogger.com/profile/17844914115422103029noreply@blogger.com0tag:blogger.com,1999:blog-7316158735602813486.post-61418581961019880742012-09-01T17:55:00.001-07:002012-09-01T17:55:43.785-07:00Sample Size Calculations for Ethics Committee Applications<br /><br />One of the key things that ethics committees look at when determining whether to accept a research proposal is the statistical justification. The ethical issue is twofold:<br /><br />1. Is an excessive number of subjects being requested while ;<br />2. Is the number of subjects sufficient to reasonably expect to get a significant result.<br /><br />Doing the statistical analysis for the application can be daunting and confusing for a prospective researcher. The online resources are often not very good either. I have found the easiest program to use for these calculations to be sigmastat. It is, however, quite expensive even for an academic license.<br /><br />Fortunately, there is a free alternative. The free software package<br />r<br /><br />can do all the power calculations for you and is available for windows, mac and even linux. It is not the most user friendly package so this guide is designed to help you use the program for your ethics commitee application.<br /><br />First decide what test you are going to use.<br /><br />For two groups with a continuous (assume normally distributed ) variable – for example the length of male vs female elephant tails – use a t test. You must also decide whether the test needs to consider whether the first group is larger or the same as the second (a one tailed test) or whether the first group could be larger, smaller or the same as the second (a two tailed test)<br />For more than one group with a continuous (assume normally distributed) variable – for example the average height of asians, africans and north americans – use an ANOVA<br />For the comparison of success or failure as a proportion of two or more groups – for example the proportion of women who conceived using a fertility treatment – use a test of proportions<br />Let’s consider these one by one<br /><br />T-test<br />work out your parameters – in this example assume that the standard deviation of the elephant tails is 20% and we hypothesize that the male ones will on average be 40% greater.<br />Fire up R and type ?power.t.test this will display the help screen<br /><br />In this case the command to get your analysis is<br /><br />power.t.test(delta=0.4,sd=0.2,sig.level=0.05,power=0.8)<br /><br />this gives you 5 subjects per group<br /><br />if you want the one tailed option (ie we are only testing whether male tails are longer male tails being significantly shorter is not a possibility)- you need<br />power.t.test(delta=0.4,sd=0.2,sig.level=0.05,power=0.8,alternative=”one.sided”)<br /><br />ANOVA test<br />This is the most tricky of the three. We first need to know our expected difference in the means and the standard deviations of the groups – is this example let’s say standard deviation again is 20% and we are looking for a difference in the means of 40%.<br /><br />The first step is to work out Cohen’s D value – this is the difference in the means divided by the standard deviation – in this case 2.0<br /><br />The next step is to work out Cohen’s effect size – this is given by<br /><br /><br /><br />so in our example we have f = 0.707<br /><br />now within r we need to load an additional module so start r and type<br /><br />require(pwr)<br /><br />then ?pwr.anova.test for the help screen<br /><br />and for our power calculation<br /><br />pwr.anova.test(k=6,f=0.707,sig.level=0.05,power=0.8)<br /><br />which gives us 6 subjects per group (you can’t use fractions of a subject)<br /><br />Test of Proportions<br />this one is fairly easy – for this example let’s say that we guess that 80% of women on a fertility treatment will conceive and of those on the placebo 10% will conceive<br />In R type<br /><br />?power.prop.test<br /><br />and for our example<br /><br />power.prop.test(p1=0.8,p2=0.1,sig.level=0.05,power=0.8)<br /><br />which gives us 7 per group<br /><br />I hope this is useful and i might revise this with some more tricky power calculations later. Any experts are welcome to comment<br />Eye Statshttp://www.blogger.com/profile/17844914115422103029noreply@blogger.com0tag:blogger.com,1999:blog-7316158735602813486.post-87544108728111790982012-09-01T17:46:00.002-07:002012-09-01T17:51:25.875-07:00Do it Yourself Climate Science with R and GHCN<br />As a scientist but yet-to-be-convinced about climate change, I was challenged to have a look at the climate data myself. After all, statistical modelling is what currently interests me and the open source software R does modelling particularly well.<br /><br />It turns out that the main surface temperature dataset upon which all of the models are based is publically available and can be analysed with R<br /><br />First you must get and install R – the program is available for all platforms from http://cran.r-project.org/<br /><br />then you need to run r and install package RghcnV3 the package written to access and transform the GHCN data<br /><br />the actual data can be downloaded from ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v3/<br /><br />monthly average, min or max data can be downloaded. This needs to be unzipped to the .dat file<br /><br />this file is then loaded into R<br /><br />> require(RghcnV3)<br />> v<-readV3Data(filename='c:\\yourdirectory\\foo.dat') #assuming your data has been renamed foo.dat and you have a windows machine - syntax is slightly different for a linux or mac<br />> stationNames<-dimnames(v) # get the station names<br /><br />Now you need the id of the station you want. This is available here http://data.giss.nasa.gov/gistemp/station_data_v2/v2.temperature.inv.txt<br /><br />for example - I tried Omeo : 50194911000<br /><br />so now which column is Omeo ?<br /><br />> which(stationNames[[1]]==”50194911000″)<br /><br />this gives column 5778<br /><br />so let’s get the data for that :<br /><br />> odata<-data.frame(v[5776,,])<br /><br />this creates a data frame of the omeo data. you can have a look at it by typing odata<br /><br />X1941 X1942 X1943 X1944 X1945 X1946 X1947 X1948 X1949 X1950 X1951 X1952<br />Jan NA NA 23.0 25.1 24.5 26.0 24.8 21.8 22.4 22.6 23.7 24.6<br />Feb NA 23.3 24.0 22.3 22.4 23.6 25.2 23.6 21.8 21.9 23.5 23.2<br />Mar NA 21.9 22.4 20.2 19.1 18.0 20.3 18.4 20.5 20.2 21.7 21.7<br />Apr NA 16.9 14.7 13.6 15.9 14.8 16.3 16.2 14.3 15.5 13.3 15.4<br />May NA 13.2 11.1 11.0 11.8 11.2 13.5 10.3 11.3 11.7 11.3 10.8<br />Jun NA 9.8 7.8 7.7 10.8 7.3 9.2 8.9 7.3 9.0 10.2 9.1<br />Jul NA 8.0 7.1 7.8 7.6 8.9 8.2 6.6 7.7 9.3 7.5 6.9<br />Aug NA 9.7 7.8 9.1 9.8 8.5 8.5 9.5 8.8 9.3 8.1 9.4<br />Sep NA 11.6 11.8 12.8 11.2 11.0 12.0 12.4 11.7 12.4 12.2 11.0<br />Oct NA 15.1 15.2 16.4 14.6 14.6 13.7 13.7 15.3 14.4 13.8 14.1<br />Nov NA 18.9 17.1 21.7 18.9 20.2 16.9 18.0 17.0 17.5 16.1 15.7<br />Dec NA 24.2 21.7 21.5 24.6 22.6 21.1 22.6 20.7 21.9 21.1 20.4<br /><br />This is a small piece of the data - in Omeo's case available mostly from 1943 to 2011.<br /><br />To get a rough and ready plot, you can stack it into a single vector<br /><br />> stackedOmeo<-stack(odata)<br /><br />and this can be plotted<br /><br />> plot(stackedOmeo$values,type=’l')<br /><br />now, there are many missing values from 1701 so let’s just plot the ones we have<br /><br />> plot(stackedOmeo$values[2900:3700],type=’l')<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-XrnuDyh0Pi0/UEKsGLmve6I/AAAAAAAAABA/raDWxTbzoz4/s1600/omeo.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="187" src="http://4.bp.blogspot.com/-XrnuDyh0Pi0/UEKsGLmve6I/AAAAAAAAABA/raDWxTbzoz4/s320/omeo.png" width="320" /></a></div><div class="separator" style="clear: both; text-align: center;"><br /></div><div class="separator" style="clear: both; text-align: left;"></div><div class="separator" style="clear: both;">and you can try and spot the warming.</div><div class="separator" style="clear: both;"><br /></div><div class="separator" style="clear: both;">R has many tools to fit a warming line : linear model, robust linear model, polynomials etc so have fun</div><br /><div class="separator" style="clear: both; text-align: left;"><br /></div><br />Eye Statshttp://www.blogger.com/profile/17844914115422103029noreply@blogger.com0