title: How to calculate if the between-subject variability is significant given repeated measurements?

text: I have a question that's been bothering me research wise. The situation is as follows: > We are measuring a quantity in subjects for which we assume that it differs > significantly between subjects, which we would like to prove. For each > subject we have 2 measurements of this quantity at two different times (time > difference between measurements is not the same for each subject). How would > I show that the variability of the quantity between subjects can not be > explained by the measurement variability? The data I have is based on voxel data, so I have for an ROI a mean and standard deviation per measurement, two measurements per subject and 50 subjects.


title: Why is the null hypothesis always a point value rather than a range in hypothesis testing?

text: This is somewhat related to another question that I asked. The question I have is, when doing hypothesis testing, when the alternative hypothesis is a range, the null hypothesis is still a point value.   As an example, when testing whether a correlation coefficient is greater than 0.5, the null hypothesis is "correlation = 0.5" instead of "correlation <= 0.5". Why is this the case? (or have I got it wrong?)


title: Probability of hitting tails

text: I have 5 variables for 1000 observations. Given the data, I need to   1. Calculate probability of hitting top 5 pct in each of the distributions separately. For example, for the value of 123.45 I need to calculate how probable is that this vale is from the top 5pct.   2. Calculate probability that a particular observation (given with 5 values for each variables) will be in the tails of each of the distributions... In the same time   3. Calculate probability that a particular observation will hit at least one 5pct tail. How would you go about this?


title: Questions about determining statistical significance in survey responses/methods tried

text: **Background:** I am studying the customer satisfaction scores from a random sampling of customers at a technology firm. I have a collection of 13 data points which consist of averages of customer satisfactions scores per year over a course of 13 years. Below is a snap shot of the years, yearly satisfactions scores, and sample size. The scoring metric ranges from 1-5, where 1 = terrible service, 5=excellent service. ![enter image description here](http://i.stack.imgur.com/vEjsm.png) **Problem:** I would like to determine if the increase in the yearly satisfaction scores are statistically significant. **Methods Tried:** I've done some googling and found the following link: "Detecting Significant Changes in Your Data", where it recommends using the T Distribution. Whenever I follow the step-by-step approach given by the author I obtain the following output, (where I average the csat scores from 2001-2003 as mu, and then average the scores from 2004-2013 as x-bar. Then I took the standard deviation of x-bar and got "s". "n" is the total number of data points since 2004, and student's t-value is calculated as t=|xbar – mu|/(s/squareroot(n)). Then p is calculated in excel as =TDIST(student t-value, degrees of freedom, tail), where those values are: 6.58, 9, and 1 respectively)). ![Output](http://i.stack.imgur.com/5yAIv.png) I obtain a value less than 5%, thus it appears that the change in csat scores is significant from 2004-2013. I have several questions about this: **Questions:**   1. Is this the correct method to use to show statistical significance in my change of csat scores? If not then what is a correct way? I should note that I have only an elementary understanding of stats, so I am limited to comfortably doing calculations in excel.   2. What is the difference between using 1-tail & 2-tails applied to my problem? Which one should I use & how do I determine that? If I use either 1-tail or 2-tails how do I determine if my % is statistically significant, would it still be the "less than 5%?" (assuming alpha=.05)   3. What is the best way to formalize my null hypothesis given this problem?   4. Why does the author (in the link) calculate p (t-distribution function in excel) using the value of "p" instead of "n"? Also, how is degrees of freedom/the tail(s) accounted for in the tdist function? Thanks for much for the help!


title: Why am I getting different intercept values in R and Java for simple linear regression?

text: I have linear regression code written in R and I have to do the same thing in Java. I used Apache Commons math library for this. I used the same data in R code and in Java code, but I got different intercept value. I could not figure out what stupid thing I have done in the code. **R Code:**               test_trait <- c( -0.48812477 , 0.33458213, -0.52754476, -0.79863471, -0.68544309, -0.12970239,  0.02355622, -0.31890850,0.34725819 , 0.08108851)     geno_A <- c(1, 0, 1, 2, 0, 0, 1, 0, 1, 0)     geno_B <- c(0, 0, 0, 1, 1, 0, 0, 0, 0, 0)      fit <- lm(test_trait ~ geno_A*geno_B)     fit      **R Output** :               Call:     lm(formula = test_trait ~ geno_A * geno_B)          Coefficients:       (Intercept)         geno_A         geno_B  geno_A:geno_B           -0.008235      -0.152979      -0.677208       0.096383       **Java Code (includes EDIT1):**               package linearregression;     import org.apache.commons.math3.stat.regression.SimpleRegression;     public class LinearRegression {         public static void main(String[] args) {                  double[][] x = {{1,0},                             {0,0},                             {1,0},                             {2,1},                             {0,1},                             {0,0},                             {1,0},                             {0,0},                             {1,0},                             {0,0}             };                  double[]y = { -0.48812477,                            0.33458213,                           -0.52754476,                           -0.79863471,                           -0.68544309,                           -0.12970239,                            0.02355622,                           -0.31890850,                            0.34725819,                            0.08108851             };             SimpleRegression regression = new SimpleRegression(true);             regression.addObservations(x,y);                  System.out.println("Intercept: \t\t"+regression.getIntercept());     // EDIT 1 -----------------------------------------------------------     System.out.println("InterceptStdErr: \t"+regression.getInterceptStdErr());     System.out.println("MeanSquareError: \t"+regression.getMeanSquareError());     System.out.println("N: \t\t\t"+regression.getN());     System.out.println("R: \t\t\t"+regression.getR());     System.out.println("RSquare: \t\t"+regression.getRSquare());     System.out.println("RegressionSumSquares: \t"+regression.getRegressionSumSquares());     System.out.println("Significance: \t\t"+regression.getSignificance());     System.out.println("Slope: \t\t\t"+regression.getSlope());     System.out.println("SlopeConfidenceInterval: "+regression.getSlopeConfidenceInterval());     System.out.println("SlopeStdErr: \t\t"+regression.getSlopeStdErr());     System.out.println("SumOfCrossProducts: \t"+regression.getSumOfCrossProducts());     System.out.println("SumSquaredErrors: \t"+regression.getSumSquaredErrors());     System.out.println("XSumSquares: \t\t"+regression.getXSumSquares());     // EDIT1 ends here --------------------------------------------------              }     }      **Java Output:**               Intercept:      -0.08732359363636362      **Java Output of EDIT1:**               Intercept:      -0.08732359363636362     InterceptStdErr:    0.17268454347538026     MeanSquareError:    0.16400973355415271     N:          10     R:          -0.3660108396736771     RSquare:        0.13396393475863017     RegressionSumSquares:   0.20296050132281976     Significance:       0.2982630977579106     Slope:          -0.21477287227272726     SlopeConfidenceInterval: 0.4452137360615129     SlopeStdErr:        0.193067188937234     SumOfCrossProducts:     -0.945000638     SumSquaredErrors:   1.3120778684332217     XSumSquares:        4.4      I will greatly appreciate your help. Thanks !


title: How to explain these p.values from model and from Anova type 2 and type 3?

text: I am trying to test the effect of each single predictor as well as their combination in predicting an outcome, and I am confusing how to explain these p values:               set.seed(123)     Y<-rnorm(100,10,3)     X1<-rbinom(100,1,0.5)     X2<-Y+runif(100,1,50)     X3<-rbinom(100,50,0.3)     fit<-glm(Y~X1+X2+X3)     #p.value(1)     round(summary(fit)$coef,4)                 Estimate Std. Error t value Pr(>|t|)     (Intercept)   8.9131     1.3292  6.7056   0.0000     X1            0.2781     0.5336  0.5211   0.6035     X2            0.0541     0.0182  2.9662   0.0038     X3           -0.0462     0.0829 -0.5577   0.5784     #p.value(2)     Anova(fit,type=2)     Analysis of Deviance Table (Type II tests)          Response: Y        LR Chisq Df Pr(>Chisq)        X1   0.2716  1   0.602278        X2   8.7983  1   0.003015 **     X3   0.3110  1   0.577077      #p.value(3)     Anova(fit,type=3)     Analysis of Deviance Table (Type III tests)          Response: Y        LR Chisq Df Pr(>Chisq)        X1   0.2716  1   0.602278        X2   8.7983  1   0.003015 **     X3   0.3110  1   0.577077       #p.value(4)     single.p.value<-NULL     for (var in c("X1","X2","X3")){         single.p.value[var]<-summary(glm(as.formula(paste("Y", var, sep="~"))))$coef[2,4]     }     single.p.value              X1          X2          X3      0.433781616 0.003440499 0.995918785       When are the `p.values` different between `car::Anova` type-II and type-III? Any body can give a clear explanation of these 4 sets of pvalues?


title: Multiple imputation using SPSS

text: I am working with a database with missing data. I have done "Roderick J. A. Little’s chi-square statistic" and knew that my data are not MCAR. However, I know don't have can I determine if there are MAR or NMAR. I was wondering if there is a test that you may suggest me. Also, how we can know that how many imputations we have to run and what factors determine the number of imputation we should apply.


title: Scale independent forecast error metric that works with changing signs

text: I am trying to analyze a quite large (~25,000 rows) dataset of cash flow forecasts. Receipts and expenses are aggregated, thus I may end up with the following data: `Forecast = 6.0e+04`   `Actual = -5.0e+04` But also with   `Forecast = 1.0e+06`   `Actual = 1.5e+06` or   `Forecast = 1.0`   `Actual = 2e+06` As you can see, the actual can differ from the forecast in order of magnitude and even signs. However, I need to find a metric for the forecast error that works for all these cases and that is scale independent. So far, I have used the Absolute Percentage Error, which works reasonably well for most of my data – but the few outliers (large forecast, small actual) render the mean absolute percentage error (MAPE) useless. I then moved to the symmetric mean absolute percentage error:   `mean(abs((act - forc)/(act + forc)))` This limits outliers as the output is between [0,1], but not if there's a sign change (=11 for the first example). Are you aware of any metrics, that limit the influence of outliers while allowing to compare across forecasting horizons and series (scale independent) and that work with changing signs? Thanks!


title: Comparing participant numbers in a systematic review

text: I am doing a systematic literature review on theories for measuring user satisfaction. I have to rate the quality of these theories, and one of the criteria should be the number of participants in the initial validation. So, we assume that each published theory was tested at least once for the initial publication (not always true, but the ones who didn't bother get a zero rating in my evaluation), and the ones who used too little participants get a rating decrease. I don't have to be very precise in this. I will be rating overall quality on a five-point ordinal scale, and I think it will be sufficient to really have only three ratings of participants, "too little", "enough" and "a whole lot", and adjust the quality rating one step up for "a whole lot" and one step down for "too little". But there is still the tripping point: when do I set the cutoff point between the three categories? Note that I am dealing with very large variations in theory and experiment design here, and with generally bad reporting. I cannot create a cutoff based on common metrics for effect size such as Cohen's d, because most of these study don't report such a metric at all. Is there any way to create rough categrization criteria across studies based on just participant number and the complexity of the theory? For the number of variables, I think that dividing the number of participants by the number of variable relations proposed by the theory will be good enough for me. (So, if study A studies the relation of one variable to satisfaction, and study B studies the relation of three variables to satisfaction, study B needs three times as many participants to be considered the same quality). But what do I do from there? I thought of making a relative comparison by using quantiles from the distributions of all studies I find, and call, say, the worst 20% of studies bad and the best 20% good. But I am not sure this will be a good strategy, as I expect a skewed distribution, and software engineering empirical experiments are known to have lower standards than those in other disciplines. I am aware that I can't do a really good comparison here, as it requires data I don't have. But is there a way of generalizing based on number of participants which is not worse than useless?


title: cluster plot: working and interpretation ?

text: Recently I have come across usage of cluster plot, which combines k-mean clustering along with PCA. The plot shows different clusters plotted using first two PCs. I have checked some of the threads (here and here) regarding the usage. I want to know, during generating a cluster plot, does the data is clustered first and then PCA is done or the reverse way (PCA followed by k-mean clustering) ? Because the second link says PCA is done followed clustering. But in the first link where an example is shown to generate a cluster plot, data is clustered first and then the cluster plot is generated. Regarding interpretation, does the plot has to be interpreted as the number of clusters generated or are there any extra points to interpret ? Please excuse me if my question is too naive. Thanks.


title: Test the randomness (uniformly distributed) on a 64 bit float random generator

text: We have a random number generator which is supposed to generate 64 bit floats, uniformly. We want to test whether it is a good uniformly random. I am not asking the general way to test it, as it was asked here http://stackoverflow.com/questions/22916519/test-the-randomness-of-a-black- box-that-outputs-random-64-bit-floats. As I understand, the first step for the test is to somehow divide the whole range into equal sizes like bins, then we keep generating and see the number of _"balls"_ falling into each bin. My question is about this first step. **How should I split the whole range of 64 bits float?** The range is between `-1.79769313486231571e+308` and `1.79769313486231571e+308` and it is huge. My idea is to utilise the **64 bits**. Can I design the equal size bins like this:   1. Every float has 64 bits, so we have 64 bins.   2. For every float generated, we read all those bits, if a bit is 1, then we increased the according bin's number by 1   3. After N samplings, the expected number of each bin should be _N/2_.    4. Then we carry on _pearson's chi-square test_ , etc.


title: How many ways to select a delegation from different sports players

text: > A school is asked to send a delegation of six pupils selected from six > badminton players, six tennis players and five squash players. No pupil > plays more than one game. The delegation is to consist of at least one, and > not more than three, players drawn from each game. Find the number of ways > in which the delegation can be selected. I use the following approach, but I failed to get the correct answer. Split the problem into seven cases:               (number of badminton players) + (number of tennis players) + (number of squash players)     1 + 2 + 3     1 + 3 + 2     2 + 1 + 3     2 + 2 + 2     2 + 3 + 1     3 + 1 + 2     3 + 2 + 1      Find the number of possible selections for each case. Add the answers, gives 10350. But the correct answer is 9450.


title: Big O notation preserved under convex functions?

text: Suppose that the random variable $X_T$ is $O_p(1)$ as $T \rightarrow \infty$. Does this imply that the random variable $\max\{0,X_T \}$ is $O_p(1)$?


title: What correlation method should I use to analyze financial timeseries?

text: I need to check if two price lists move togheter. I three methods:   1. Pearson   2. Spearman   3. Kendall Could someone explain me why I should use a method instead of other to check if those price lists move in the same direction? I know that Pearson is often used for this kind of analysis... Is that correct? Thank you


title: How to formally test for a "break" in a normal (or other) distribution

text: It frequently comes up in social science that variables that _should_ be distributed in some way, say normally, end up having a discontinuity in their distribution around certain points. For instance, if there are specific cutoffs such as "passing/failing" and if these measures are subject to distortion, there may be a discontinuity at that point. One prominent example (cited below) comes student standardized test scores are normally distributed basically everywhere except at 60% where there is very little mass from 50-60% and an excessive of mass around 60-65%. This occurs in cases where teachers grade their own students exams. The authors investigate whether teachers are really helping students pass exams. The most convincing evidence without a doubt comes from showing the graphs of a bell curve with a big discontinuity around different cutoffs for different tests. However, how would you go about developing a statistical test? They tried interpolation and then comparing fraction above or below and also a t-test on the fraction 5 points above and below the cutoff. While sensible, these are ad-hoc. Can anyone think of anything better? Link: **Rules and Discretion in the Evaluation of Students and Schools: The Case of the New York Regents Examinations** http://www.econ.berkeley.edu/~jmccrary/nys_regents_djmr_feb_23_2011.pdf ![Distribution of test scores, manipulable ones in black, note the sharp drop in density below the cutoff and corresponding rise above](http://i.stack.imgur.com/QYdGN.jpg)


title: How do I mathematically prove that k-means clustering converges to minimum squared error?

text: I am using k-means clustering to analyze and obtain patterns in traffic data. This well-known algorithm performs 2 steps per iteration.   1. Assign each object to a cluster closest to it, based on the euclidean distance between the object and a cluster's center (obtained at a previous iteration).   2. Update the centers of the clusters based on the new members of clusters. I understand that k-means clustering overall aims to minimize the sum of the square errors. However, how can we mathematically prove that the 2 steps eventually converges to this goal? It sounds like a complicated optimization problem involving Lagrange multipliers, but maybe there is a more straightforward way to prove it?


title: Algorithms to eliminate or display outliers in a dataset

text: I am trying to visualize a large multi-dimensional dataset. Each dimension has a different range. Values in one column may range between 0-100 while values in another could range from a few hundred thousand to a few hundred million. Therefore it is really hard to show a graph with a reasonable scale. I would like to visualize them using a parallel coordinate chart or series of point charts. I have come up with three different approaches; any outlier outside of the standard deviation will be removed, any outlier outside of the standard deviation will show up at the edge of the standard deviation, or I will adjust the scale of each dimension so that all are scaled to be seen at once. Does anyone have any recommendations regarding which method would retain the integrity of the data? Also, does anyone have good examples of outlier removal algorithms?


title: How well do analyst ratings predict stock perfromance?

text: Consider a huge dataset of which the following is a characteristic sample:               stock.id date analyst1 analyst2 analyst3 return            1 Jan         2       -2        1     10            2 Feb        -1        2        0    -10            3 Mar         0        2        0      3            4 Apr         0        0        0      2            5 May         1        0        2     -1            6 Jun         1        0        2     -3          ...      The "analyst" columns are ratings issued by individual stock analysts and the "return" column is the subsequent return on rated stock. Here rating 2 means a strong buy and -2 means a strong sell. I am looking for a statistical test or techniques that can tell me how well the analysts performed. Ideally, I want to come up with scores to indicate which analysts have rated stocks correctly.


title: Relations between distribution free statistics and nonparametric statistics?

text: Is a distribution free statistic necessarily a nonparametric statistic? Is there a distribution-free statistic that is not nonparametric? Is a nonparametric statistic necessarily a distribution free statistic? Is there a nonparametric statistic that is not distribution-free? Are distribution free statistics and nonparametric statistics (both viewed as measurable mappings of sample points) the same concept? How is a nonparametric statistic defined? Is a nonparametric statistic defined as a statistic whose distribution doesn't depend on the distribution of the sample? Thanks and regards!


title: Interpreting weird box plot with reversed whiskers

text: Although I'm fairly new to box plots, I thought I had got the hang of them, until I came upon this one today. I don't know what to make of that bottom "whisker" drawn inside the box. ![enter image description here](http://i.stack.imgur.com/2JQ3O.png) This population is composed by only four values: 16.5, 17.14, 13.5, 16.75 Granted, the small size of this population is probably not ideal for this kind of diagram. But this is just one among several ones that I'm plotting, and the others make a bit more sense. My question is two-fold:   1. Is this even a valid representation or is the software that I'm using to draw it misbehaving? (I'm using python and matplotlib)   2. If it is valid, how do I go about interpreting it?  * * * **Edit:** Including a figure created with a `whis=3` parameter (see answers below): ![enter image description here](http://i.stack.imgur.com/A9LC8.png)


title: Derivation of uncertainty propogation?

text: For a function $$ y =f(x), x=\left(x_1, x_2, ..., x_N\right)$$ the law of propagation of uncertainty, see GUM sect 5$^{[1]}$ (pdf), is generally given as $$ u_y^2 = \sum_{i=1}^N \left(\frac{\partial f}{\partial x_i}\right)^2 u_{x_i}^2 $$ where $u_{x_i}$ is the uncertainty associated with each $x_i$. How do we derive this propagation law? I understand that it can be achieved by Taylor expanding the function, but I don't see why this is a sensible thing to do or what is the physical meaning of it. (I've asked basically the same question at Signal Processing here but haven't got any useful answers. Hopefully you guys are more suited to this question). [1] _Evaluation of measurement data – Guide to the expression of uncertainty in measurement_   JCGM 100:2008 (html) [Sec 5 here]


title: why is orthogonality represented geometrically in DOE and is fractional factorial designs are orthogonal?

text: I've noticed in many places that it is common do describe experiments in DOE and orthogonality as a geometric figure, why is that? Also is it connected to orthogonality?


title: Needing help in order to infer the statistical hypothesis tests performed in an old paper

text: I am in great need of help in order to infer the statistical hypothesis tests performed in an old paper. However I need to make some reasonable guesses only from the abstract of the paper since the original text is in chinese which i can not understand (after extensive search i was also unable to find the original text even in chinese) The title of the paper is: **“Detection of sister chromatic exchange in workers exposed to coal tar pitch and to coke oven volatiles"** The abstract of a paper (which was published in 1998 is the following): In order to know the changes of genetic toxicological effects on workers occupationally exposed to polycyclic aromatic hydrocarbons (PAHs), sister chromatic exchange(SCE) was detected by the methods of peripheral lymphocyte culture in 23 workers exposed to coal tar pitch (CTP) and in 19 workers exposed to coke oven volatiles (COV) and 12 normal controls. The results suggested that the SCE in occupational workers was significantly higher than that in controls (11.31 vs 6.37, P < 0.001). The SCE in workers exposed to CTP and to COV was higher than that of control (10.27 and 12.58 vs 6.37) respectively. In workers exposed to CTP and COV, there were no differences of SCE for smokers and nonsmokers (P > 0.05). It is indicated that CTP and COV caused strong genetic toxicity and injury to chromosome. In your opinion how do you think that the above reasarch was organized For example: a) What types of statistical hypothesis testing was performed by the reasearchers b) What kind of data was collected and used for each statistical hypothesis test c) What methods were employed for the each hypothesis test ?


title: Interactions between random effects

text: I'm considering a mixed-effects model to try to understand factors that influence the number of ticks sampled on wild rodents. My data is nested so that I have one tick count per rodent, multiple rodents per site and multiple sites per year (sites are repeated across years but not every site is present in all years). So far, without fixed effects, my model looks something like this:               glmer(Ticks~(1|Year)+(1|Site), family=poisson, data=tickdata)      I understand that this will account for random variation between sites and between years. My main concern is that I think the effect of specific sites will change between years, e.g. some sites may recieve more rainfall than others, but the identity of the sites with the highest rainfall differs between years. So do I need some sort of interaction term? Maybe something like this:               glmer(Ticks~(1|Year)+(1|Site)+(1|Year:Site), family=poisson, data=tickdata)       Thank you in advance.


title: Visualize effects of hurdles on user numbers

text: I have the following data:               1000 users were invited      990 emails delivered successfully      200 visited the system       60 tried to create an account        58 tried to create an account with a valid email       50 successfully created an account (i.e. chose a secure password)       40 participated by voting       10 posted something      While I could use a bar graph, the proportions of 60:58:50 would likely be lost. If I'd use a logarithmic scale, the proportions would be distorted and the feeling of numbers would be lost. Is there a good way to represent both the (significant) loss in user attention due to the number of hurdles as well as the effect of each hurdle?


title: Testing the difference of (some) quantile-Q between groups?

text: For some Y variable, that is divided into 3 groups (X), I wish to compare the groups and for the hypothesis that the 90% quantile is the same between all three groups. What tests can I use? One option I can think of is using quantile regression, are there other alternatives/approuches? I imagine that if I had wanted to compare the median, I could have used the kruskal wallis test (though it is based on ranks, but if I remember correctly, it would give the same results when the residual distribution is symmetric) Thanks.


title: Why is the sampling distribution of normal distributed variable automatically also normal distributed

text: I am currently reading about the standard error. I know about central limit theorem, but I don't understand why, if my variable is normally distributed in the population, the sampling distribution (the distribution of the sample means) is also a normal distribution, no matter of sample size etc.


title: Instability in glmer model at zero random effect variance

text: I am getting very different results from glmer (using version 0.999375.37 of lme4) when running in different environments (I can't pinpoint the source of discrepancy, but let's assume some floating point precision difference in the input data). Here is the invocation:               glmer(data=data, n ~ 1 + offset(log(mu)) + (1 | id),       family=poisson(link=log), weight=w,verbose=T)      Sometimes the output is this:                 0:     532.77427: 0.610172 -0.0811043       1:     525.91180: 0.531918 -0.102716       2:     511.77820: 0.373325 -0.137516       3:     495.51600:  0.00000 -0.127691       4:     478.25662:  0.00000 -0.100781       5:     474.47031:  0.00000 -0.0806405       6:     474.46819:  0.00000 -0.0811088       7:     474.46819:  0.00000 -0.0811044       8:     474.46819:  0.00000 -0.0811043       9:     474.46819:  0.00000 -0.0811043     Generalized linear mixed model fit by the Laplace approximation      Formula: n ~ 1 + offset(log(mu)) + (1 | id)         Data: data               AIC           BIC        logLik    deviance      478.4681939 488.556259695 -237.23409695 474.4681939     Random effects:      Groups Name        Variance Std.Dev.      id     (Intercept)  0        0           Number of obs: 1146, groups: id, 160          Fixed effects:                         Estimate       Std. Error  z value   Pr(>|z|)         (Intercept) -0.0811043131501  0.0100762538900 -8.04905 8.3436e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      and sometimes it's this:                 0:     532.77427: 0.610172 -0.0811043       1:     525.91180: 0.531918 -0.102716       2:     511.77820: 0.373325 -0.137516       3:     495.51540:  0.00000 -0.127690       4:     478.25643:  0.00000 -0.100780       5:     474.47031:  0.00000 -0.0806405       6:     474.46819:  0.00000 -0.0811088       7:     474.46819:  0.00000 -0.0811044       8:     474.46819:  0.00000 -0.0811043       9:     474.46819:  0.00000 -0.0811043      10:     474.46819:  0.00000 -0.0811043      11:     474.46819: 1.29719e-05 0.568362     Generalized linear mixed model fit by the Laplace approximation      Formula: n ~ 1 + offset(log(mu)) + (1 | id)         Data: data       AIC  BIC logLik deviance      5700 5710  -2848     5696     Random effects:      Groups      Name        Variance   Std.Dev.        id          (Intercept) 1.6827e-10 1.2972e-05     Number of obs: 1146, groups: id, 160          Fixed effects:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) 0.568362   0.007282   78.05   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Note that fitting the model without the random effect gives:               Call:  glm(formula = n ~ 1 + offset(log(mu)), family = poisson(link = log),          data = data, weights = w)          Coefficients:          (Intercept)       -0.0811043131378        First question: Why is the result unstable? IIUC the model is of the form n ~ Pois(mu * exp(a + b_id)) where the b_id are drawn from a normal distribution N(0, theta). In the verbose output the third column is a, and the second column is related to theta. As theta approaches zero, the b_id's must go to zero, so I would expect a to converge to the maximum likelihood estimate for the model without random effects. This is exactly what we see, except for iteration #11 in the second outcome. Second question: How can I automatically detect and avoid this? Is it sufficient to fall back to the simple model (without random effect) when theta is below some threshold? Note that I'm interested in the estimated values for each level, not the variability (the opposite of the situation discussed in this FAQ). I also want some penalty on differences between the values for different id's, so I don't want to model id as a fixed effect.


title: What is name of this distribution and can we have 68–95–99.7 rule for it?

text: I have a distribution like this: ![enter image description here](http://i.stack.imgur.com/wHaux.jpg)   * What is name of this distribution?   * As you know we have `68–95–99.7` rule in `normal distribution`. Can we have something like this in this distribution? Thanks.


title: Which model should I use to fit my data ? ordinal and non-ordinal, not normal and not homoscedastic

text: Here is the kind of data I have: I have two predictor variables: 1) discrete non-ordinal --> c('a','b','c') 2) discrete ordinal --> c(10,100,200,500) Response variable: Proportion of TRUE over a list of TRUE/FALSE. If the list is of length 3, my variable can take only 4 values. But not all my values come from the same list's length. And moreover the lists are globally long ! So it is discrete proportions but can take more than 100 values. Here is an example (resp is a subset of my data):               pred_1 = rep(c(10,20,50,100),30)     pred_2 = rep(c('a','b','c'),40)          resp = c(0.08666667, 0.04000000, 0.13333333, 0.04666667, 0.50000000, 0.04000000, 0.02666667, 0.24666667, 0.15333333, 0.04000000, 0.06666667, 0.06666667, 0.03333333,         0.04000000, 0.26000000, 0.04000000, 0.04000000, 1.00000000, 0.28666667, 0.03333333, 0.06666667, 0.15333333, 0.06666667, 0.28000000, 0.35333333, 0.06000000,         0.06000000, 0.05333333, 0.96666667, 0.06666667, 0.03333333, 0.22000000, 0.04666667, 0.04666667, 0.05333333, 0.05333333, 0.05333333, 0.08000000, 0.48666667,         0.08666667, 0.02666667, 0.21333333, 0.45333333, 0.04666667, 0.36000000, 0.06666667, 0.04000000, 0.06000000, 0.07333333, 0.06000000, 0.04000000, 0.04666667,         0.30000000, 0.08666667, 0.07333333, 0.06666667, 0.29333333, 0.36000000, 0.17333333, 0.04000000, 0.09333333, 0.11333333, 0.03333333, 0.08000000, 0.27333333,         0.08666667, 0.03333333, 0.04000000, 0.02666667, 0.07333333, 0.07333333, 0.02000000, 0.02666667, 0.08000000, 0.07333333, 0.02666667, 0.06666667, 0.07333333,         0.95333333, 0.05333333, 0.04000000, 0.11333333, 0.04000000, 0.07333333, 0.06666667, 0.05333333, 0.04000000, 0.04000000, 0.06000000, 0.12666667, 0.04666667,         0.04000000, 0.21333333, 0.05333333, 0.97333333, 0.11333333, 0.02666667, 0.04000000, 0.03333333, 0.37333333, 0.25333333, 0.06000000, 0.06000000, 0.06000000,         0.04666667, 0.26666667, 0.98000000, 0.02000000, 0.26000000, 0.06000000, 0.05333333, 0.28000000, 0.99333333, 0.04666667, 0.02666667, 0.04000000, 0.12666667,         0.04666667, 0.18000000, 0.03333333)       my response variable is not at all normally distributed (kolmogorov-smirnow and shapiro test + visual checking with qqplot()) nor is the residuals of a linear model (lm()). Moreover the common assumption of homoscedasticity is not respected neither. I've always asked a similar question but not as much accurate here. Peter Flom has suggested that I use a ordinal logistic regression (polr()) but I might not have given him enough information (he did not know the number of levels for example). What do you think ? Which model would you suggest me ? Can I make a polr() with that much levels ? When I do it I actually get this: Error message: "Initial value "vmin" in not finite" Notification message: "glm.fit: fitted probabilities numerically 0 or 1 occurred " I'm struggling on this problem for quite a long time. All your contributions are more than welcome ! Thanks a lot !


title: Post hoc tests for robust mixed design ANOVA using R

text: Is it possible to compute the function `mcp2atm` with unequal sample size? I ran the robust Mixed Design ANOVA by using `tsplit` and want to do post hoc tests now.


title: How can a multiclass perceptron work?

text: I don't have any background in math, but I understand how the simple Perceptron works and I think I grasp the concept of a hyperplane (I imagine it geometrically as a plane in 3D space which seperates two point clouds, just as a line separates two point clouds in 2D space). But I don't understand how one plane or one line could separate three different point clouds in 3D space or in 2D space, respectively – this is geometrically not possible, is it? I tried to understand the corresponding section in the Wikipedia article, but already failed miserably at the sentence “Here, the input x and the output y are drawn from arbitrary sets”. Could somebody explain the multiclass perceptron to me and how it goes with the idea of the hyperplane, or maybe point me to a not-so-mathematical explanation?


title: MCMC autocorrelation convergence diagnostic

text: I use MCMC (Metropolis-Hastings) to sample posterior distributions of three parameters using a nonlinear least-squares objective function to calculate the likelihood of a parameter sets. I want posteriors for about 500 different conditions, each time I calculate a posterior with MCMC. Since it is hard to check all 500 runs by hand for convergence, I would like to have a simple rule based on autocorrelation which indicates for each run if the chain converged or not, checking all conditions individually would be very cumbersome. Each chain has length 40000 (burn-in already removed) and I calculated autocorrelations for lags 1 to 10000. I then calculated the autocorrelation time as the time constant for the autocorrelation series (the first lag when the autocorrelation drops below `1/exp(1)`). I get autocorrelation times between 16 and about 7000. I suppose that with an autocorrelation time of 16 I can be sure that the Markov chain converged, but is there maybe some "rule of thumb" which tells me, given a chain of length 40000, which autocorrelation time is still "acceptable"? I searched and looked for instance here, however, the approaches I found require a visual investigation.


title: Measuring forecast accuracy

text: We're forecasting sales data for one of our clients on a weekly basis. Sales is forecasted for each organizational unit. The sales data is forecasted via different algorithms and/or algorithm parameters. So, for every organizational unit, about 10 different forecasts are generated. Now I want to find the best forecast/algorithm to be chosen as the weekly forecast based on the accuracy of each algorithm of the past 9 weeks. So for each algorithm there is data like:               Org Unit     Algorithm      Date          Actual       Forecast     -----------------------------------------------------------------     OU1          A1             2013-07-25    100          110       OU1          A1             2013-07-24    100          120       OU1          A1             2013-07-23    130          130       OU1          A1             2013-07-22    140          170       OU1          A1             2013-07-21    110          130       OU1          A1             ...           ...          ...     OU1          A2             2013-07-25    100          102     OU1          A2             2013-07-24    100          108      OU1          A2             2013-07-23    130          120       OU1          A2             2013-07-22    140          122       OU1          A2             2013-07-21    110          130       OU1          A2             ...           ...          ...      Based on the sample data shown above I need to decide whether to choose A1 or A2. What would you suggest? Standard deviation on the mean absolute error plus average of MAE? Theil's U? P.S. It does not matter whether the forecast is below or above the actual, only the absolute deviation needs to be considered.


title: Multivariate analysis for factors affecting survival

text: I am doing an audit looking at survival for a particular cancer. I have done a univariate analyses and Kaplan-Meier curves using Prism Graphpad looking at prognostic factors which may affect survival. I now want to enter the significant factors into a multivariate analysis. I shall be using SPSS version 21. My question is that one of my factors, histology showed a trend for significance and I would like to include that. But there are 5 options for different histologies, and there is no reference category for that. I have managed to do some of my other multivariate Cox regression on SPSS, via Analyze, Survival, Cox regression. Is it possible to do it the same way? Is this a matter of using a specific coding of categorical variables?


title: How to Analyse and predict

text: I'm having a set of data's, say               Question,     Question's Main Category,Main Category followers,     Question's Related Categories,Each Related Categories followers,     Whether it got answer in 24hrs [Yes/No]       Now, I need to analyse a database with thousands of sets of data's like given above and predict whether the new given question(with all details like main category,followers...etc) will get answer in 24hrs or not ? I have no idea which algorithm or statistical analysis will be more suitable to solve this kind of problem.


title: Multiple testing, FDR and power

text: Suppose that I've a data with a total of 30,000 observations. Each observation has a total of 6 values, 3 per condition. I'm interested in, let's say, testing if the means of these two conditions for each of these observations are different or not. For this, we typically apply a statistical test of choice to each of these 30K observations. In order to control the false positive rate due to multiple testing, we'd control for FDR at say 5% using BH method. Example of an observation would be:               # dummy data     Condition 1: 40, 55, 48     Condition 2: 129, 77, 181      Also, let's assume I obtained a total of 1000 "significant" results. If I _filter_ the data somehow, say, _remove all those observations where the maximum value across all 6 values is < 5_ and end up with a filtered set of 20K observations, then I find that the number of significant results is a bit higher, say 1250. If I continue this approach of filtering data with higher and higher thresholds (< 10, < 20 ... preferably using quantile = 0.1, 0.2, 0.3 etc..), I find that the number of significant results at 5% FDR keeps increasing up to a certain point and with _too stringent filtering_ starts to reduce again. It's fairly obvious that in controlling the FP in a multiple testing setup one compromises for the statistical power. My question is, are there methods that could somehow compute a filtering criteria which maximises the power (to detect after FDR correction)? If not, is it statistically sound approach to try and filter data with more than 1 value and decide on the one that maximises significance?


title: The result of absolute change and relative change in t test is different

text: My task is to evaluate the effectiveness of an intervention. In my study, there are one control group and one intervention group. There are two time points, baseline and endpoint. And I use the score of some questionnaire as outcome measures. My strategy to show/test the effectiveness is to compare the change between baseline and endpoint between these two groups by independent t test for continuous measures. I calculated both the point difference $(\text{endpoint}-\text{baseline})$ and relative difference $((\text{endpoint}-\text{baseline})/\text{baseline})$. After running t test in SPSS, I got two different result that the $p$-value of the $t$-test on point difference was higher than $0.05$ and $p$-value of the $t$-test on relative difference is lower than $0.05$. These results were contradicting because I set the alpha as $0.05$. I want to ask: **which difference should I use?**


title: Combining two sets of principal components in a regression

text: My data set consists of a large number of time series of two frequencies: daily and monthly. I am interested in doing a predictive regression on a monthly and weekly basis. I am thinking of dividing my data set in two (daily and monthly), and then do principal component analysis on each subset individually and use the two separate sets of principal components in my predictive regression. Are there any contraindications in combining two separate sets of principal components (one from the daily and one from the monthly data) in a regression?


title: Split testing significance with continuous measures

text: What technique would you use to understand the statical significance of a test where the target outcome is a continuous measure? Say there are three tests Control, A and B each with 50k samples and they generate the following revenues of 1054, 750 and 450. I've done this with discrete measures like 20 people out of 1000 converted but never continuous measures.


If you flip a fair coin indefinitely, how likely is it that, at some point, you will have evidence that the coin is not fair? What I mean is this: Take a fair coin and select a $p$-value. If we were to flip the coin $100$ times, there would be some threshold proportion $Q$ such that if the propotion of heads (or tails) were larger than $Q$, we would have evidence (given our selected $p$-value) that the coin was not fair. If we were to flip the coin $1000$ times, the threshold proportion would diminish. As the number of flips $N$ increases without bound, the $Q$ approaches $.5$. If we flip the fair coin indefinitely, what is the probability that there will be some $N$, such that, if we had stopped flipping at $N$, we would have ended up with evidence against the coin's being fair (in terms of a $p$-value)?


title: Likelihood analysis for exponential distribution

text: Assume a collection of independent exponential random variables $y_{1}, \ldots, y_{n}$ with means $\mu_{1}, \ldots, \mu_{n}$; where $\mu_{i} = \beta_{0}+\beta_{1}x_{i}$. How can I find the profile likelihood of $\beta_{1}$? From what I gather, the profile likelihood is a technique for eliminating so- called "nuisance parameters" from a given inference. This is typically achieved by maximising the likelihood with respect to the nuisance parameters whilst holding the parameters of interest fixed. For the two-parameter model described above, because we are only interested in $\beta_{1}$, we can consider $\beta_{0}$ as a nuisance parameter. **ATTEMPT:** Firstly, given that $\mu_{1}, \ldots, \mu_{n}$ are the means for $y_{1}, \ldots, y_{n}$ exponential variables, and, in general notation, the mean of the exponential distribution is given by $\lambda^{-1}$, then each random variable must follow an $\text{Exp}((\lambda^{-1})^{-1}) \implies \text{Exp}(\mu^{-1})$; in other words: For $i=1, \ldots, n$: $$ \mu_{i}^{-1}\text{exp}\left\{-\mu_{i}^{-1}y_{i}\right\} $$ Now, to determine the form of the likelihood, we form $p(x_{1}) \times p(x_{2}) \times \ldots \times p(x_{n})$, hence, $$ L(\mu) = \prod_{i=1}^{n} \left[\mu_{i}^{-1}\text{exp}\left\{-\mu_{i}^{-1}y_{i}\right\}\right] $$ Simplifying we get, $$ = \left[\prod_{i=1}^{n}\mu_{i}^{-1}\right]\text{exp}\left\{-\sum_{i=1}^{n}\frac{y_{i}}{\mu_{i}}\right\} $$ Because $\text{log}\left(\prod_{i=1}^{n}\mu_{i}^{-1}\right) = \sum_{i=1}^{n}\text{log}\left(\mu_{i}^{-1}\right)$, we get: $$ \text{log}(L(\theta)) = \sum_{i=1}^{n}\text{log}\left(\mu_{i}^{-1}\right)-\sum_{i=1}^{n}\frac{y_{i}}{\mu_{i}} $$ Now, subbing in for $\mu_{i} = \beta_{0}+\beta_{1}x_{i}$, gives: $$ = \sum_{i=1}^{n}\text{log}(\beta_{0}+\beta_{1}x_{i})^{-1}-\sum_{i=1}^{n}\frac{y_{i}}{\beta_{0}+\beta_{1}x_{i}} $$ Now that we have the log-likelihood, can we just find the MLE for $\beta_{0}$, substitute the acquired MLE for $\beta_{0}$ back into the log-likelihood and simplify? Will this achieve the profile likelihood for $\beta_{1}$? In other words, where $\hat{\beta_{0}}$ is the MLE for $\beta_{0}$, will the profile likelihood be given by the following? $$ L_{p}(\beta_{1}) = \sum_{i=1}^{n}\text{log}(\hat{\beta}_{0}+\beta_{1}x_{i})^{-1}-\sum_{i=1}^{n}\frac{y_{i}}{\hat{\beta}_{0}+\beta_{1}x_{i}} $$ **Note the "hat" notation for $\beta_{0}$in the above.**


title: Model matrix for quadratic regression with variable number of covariates in R

text: I want to create a model matrix for a quadratic regression in R. I have a matrix $X$, where each column is a covariate. Currently, I'm using the code:               Model = model.matrix(~ 1 + X[,1] + X[,2] + X[,3] + I(X[,1]^2/2) +      I(X[,2]^2/2) + I(X[,3]^2/2)  + I(X[,1]*X[,2]) + I(X[,1]*X[,3])+I(X[,2]*X[,3]))      Then if I'm dealing with a problem with 4 or more covariates I have to add the new terms by hand. I tried to create the model matrix my self using:               Model = cbind(1,X,X^2/2,mixed_terms)      I don't remember exactly how I was calculating the mixed terms (I think I used a permutation of the columns of $X$), but in any case my code was much slower the model.matrix, while I need the model matrix to be create fast because it goes inside a loop in simulation routine. Thanks!


title: In a time series study, how many time points are enough for unit root test

text: I'm keeping some fish for toxicity test and will analyze specific parameter after exposing toxicants. I am keeping the fish in normal condition and check the parameter daily as "baseline". I plan to run unit root test to confirm whether there is change after toxicant exposure. However, I'm not quite sure how long the "baseline" data I need to obtain will considered "enough" for the upcoming analysis. Anyone may provide me hints? Thank you.


title: Reporting predicted outcome of mixed effects models

text: I am analyzing some group-randomized longitudinal trial data with three time points using mixed effects model. Some people missed 1-2 time points, and I didn't exclude any cases with missing dependent variable prior to running the analysis. When I was asked by a colleague to report the predicted value, I hit into a problem: Let's say the simple form of the regression model is: $y = \beta_0 + \beta_1 x1 + \beta_2 x2$ After estimating the regression coefficient, we can then find out the predicted y, $\hat{y}$. My question is, when calculating this predicted y, should I include cases that have a full set of predictors, but missing on the outcome? Is there any conventional rule to follow? Thank you very much for your advice in advance.


title: How is the Df computed in a mixed model?

text: The following is the output for a mixed model example. The only difference between fm1 and fm2 is the random factor "URBAN", why the df for fm2 is 5 but not 4? Any help would be great.               > anova(fm.null, fm1, fm2)     Data: example     Models:     fm.null: SCIENCE ~ 1 + (1 | GROUP)     fm1: SCIENCE ~ URBAN + (1 | GROUP)     fm2: SCIENCE ~ URBAN + (URBAN | GROUP)             Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)     fm.null 2 644.17 650.32 -320.08     fm1     3 505.37 514.60 -249.69 140.79 1 < 2.2e-16 ***     fm2     5 423.22 438.60 -206.61 86.15 2 < 2.2e-16 ***     ---     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


title: Why is the Likelihood function NOT a case of the inverse fallacy?

text: This may be a trivial question, but as a research psychologist I do not have a robust statistics background to answer it. It appears to me that the likelihood function--$L(\theta | \text{data}) = P(\text{data} | \theta)$--is committing the inverse fallacy which is exactly what using Bayes' theorem avoids. I'm sure that the logic behind the likelihood is sound, but I can't see why this is NOT a case of incorrectly equating two different conditional probabilities (i.e., the inverse fallacy).


title: Using Rasch model to explain relationships between a set of dependent and independent variables

text: My research study is in development economics. My data consist of more than one independent variables (continuous and categorical) as well as more than one dependent variables (categorical 5-point Likert scale). I want to find out their relationships, therefore regression (ordinal logistic) is the best choice. But the problem is to reduce the dependent variable from many (6 variables) to one. I explored many options such as categorical principal component analysis (CATPCA), principal component analysis (PCA), composite variable (COMP). But using PCA and COMP means considering the data as continuous instead of categorical. However, all lead to a reduced dependent variable (object dimension/factors/composite) that is in continuous data form, for which I have to use linear regression. If I choose ordinal logistic regression then I have to run it on each dependent variable separately, leading to many models instead of one model for dependent variable. I just came across the Rasch model. I am not much familiar with it. It seems to fit my data. My questions are:   1. Could I use Rasch Model?   2. If yes then which one, the rating scale or partial credit in rasch model OR Graded response model?   3. Could I able to get results from it on which I can run regression (linear/ordinal?)?   4. Which software is best for it means easy to use and free to download?


title: Is listwise deletion / complete case analysis biased if data are not missing completely at random?

text: In the comments to the answer to my question I stated _"Many rows have only 1 missing variable, so to exclude the row think leads to bias (they are not MCAR)"_ and in reply I was told _"You're wrong, see Rubin's Statistical Analysis with Missing Data 2nd ed. CC is unbiased with MAR data,"_ I don't have Rubin and Little's book but I had been able to borrow it a few months ago and I am sure I learned that complete case analysis is biased unless the data are missing completely at random. Now I am terribly confused. Can anyone explain / clarify my misunderstanding ?


title: Matching curves

text: I've a large number of curves over a finite interval. Given a new curve C, I want to find a set of curves in the database s.t the area between a matched curve and C is less than a certain threshold. Is this a standard problem? If so, what is it called? What algorithm can I use to solve this (besides the trivial algorithm of iterating through all curves)?


title: Confidence intervals on differences in choices in a GEE framework: methods and alternatives?

text: Say we have 5 items, and people are asked which item they like. Multiple answers are possible, but also no answer is possible. The people are categorized according to factors like gender, age, and so on. One possible approach to analyze the differences between genders, age groups and so on is making use of the Generalized Estimating Equations. I thus construct a dataset with a binary variable indicating whether the item was liked or not, and as predictor variables I have the items, the person id, the age,... i.e. :               Liked   Item   Person ...     0       1      1     1       2      1     0       3      1     0       4      1     1       5      1     1       1      2     ...      Then I apply a model with following form : $$Liked = Item + Gender + Item*Gender + Age + Item*Age + ... $$ with Person as random factor (called id in some applications) and a logit link function. Now I like to give confidence intervals on the conditional fractions, i.e. the confidence interval of the fractions males and females that liked a particular item, corrected for age differences and the likes. I know I could use the estimated coefficients to get the results I want, but I'm a bit lost in how to do that. I can calculate the estimated odds, but I'm not sure on the standard error (SE) on those odds based on the SE of the coefficients. I'm not sure on how to deal with the random component of the variance for example. So: 1) Any pointers on how to calculate that SE from the SE of the coefficients? 2) Any alternatives for an approach? I've been thinking about mixed models, but a colleague directed me to GEE as more appropriate for these data. Your ideas? * * * Edit : for practical pointers, I'm using geepack in R for this. I tried `effect()`, but the option `se.fit=T` is not implemented. In any case, that would give the SE for every observation, which is not what I'm interested in.


title: For PCA prove that covariance matrix is diagonal?

text: For $z=W^\mathbf{\intercal}(x-m)$ where $k$ columns of $W$ are the top $k$ principal components of a given dataset $x$ and $m$ is the mean, how can I prove that $z$ is diagonal? $x$ is the original dataset (dxN dimensional matrix) and $m$ is the mean of $x$(1xN dimensional matrix) and $k$ is a number less than d


title: Cox Baseline Hazard

text: Let's say I have a "kidney catheter" data set. I'm trying to model a survival curve using a Cox model. If I consider a Cox model: $$h(t,Z) = h_0 exp(b'Z),$$ I need the estimate of the baseline hazard. By using the built-in `survival` package R function `basehaz()`, I can easily do it like this:               library(survival)          data(kidney)     fit <- coxph(Surv(time, status) ~ age , kidney)     basehaz(fit)      But if I want to write a step by step function of the baseline hazard for a given estimate of parameter `b` how can I proceed? I tried:               bhaz <- function(beta, time, status, x) {              data <- data.frame(time,status,x)         data <- data[order(data$time), ]         dt <- data$time         k <- length(dt)         risk <- exp(data.matrix(data[,-c(1:2)]) %*% beta)         h <- rep(0,k)              for(i in 1:k) {             h[i] <- data$status[data$time==dt[i]]/sum(risk[data$time >= dt[i]])                   }              return(data.frame(h, dt))     }          h0 <- bhaz(fit$coef, kidney$time, kidney$status, kidney$age)      But this does not give the same result as `basehaz(fit)`. What is the problem?


title: How do Bayesian GLMs with noninformative priors on the coefficients compare to estimation using MLE's?

text: I'm curious as to how putting noninformative priors on the regression coefficients in a GLM compares to maximum likelihood estimation, _in frequentist terms_. A dispersion parameter $\phi$ is of course a source of trickiness, so perhaps we assume that these are estimated using empirical Bayes. In the simplest situation, the response $Y_i, i = 1, ..., n$ might be Gaussian distributed with mean $$ E(Y_i | \beta) = x_i ^T \beta $$ have common variance $\sigma^2$, and be independent (given $\beta$). Then the Bayes estimator of $\beta$ is the same as the MLE, and for $\sigma^2$ estimation via empirical Bayes results in $\hat \sigma^2 = \frac{\sum(Y_i - \hat Y_i)^2}{n - p}$ where $p$ is the number of $\beta$'s, which makes an appropriate correction to the denominator that the MLE does not make. My general sense is that it is a bit easier to get away with flat priors on mean parameters, and that integrating out the mean parameters is generally what you want to be doing to estimate the dispersion parameter $\phi$ (short of putting an informative or "weakly informative" (in Gelman's sense) prior on it). That is to say, we would rather integrate out $\beta$ than max- marginalize (aka profile) it out. It isn't obvious to me immediately that flat priors on regression coefficients lead to proper posteriors in GLMs; part of me feels that they usually ought to, but part of me also feels that maybe they run into the same sort of issues as MLE's; for example, maybe flat priors break down if the data are perfectly seperated in a logistic regression. So, general insight or guidance on this is appreciated, as well as references.


title: Estimator for an incidence rate

text: While going through a statistics course for medicine students, I ran accross a problem related to incidence rates. The context of the problem is a chapter about the Poisson distribution. In the problem, 2300 smokers are followed over a 1 year span during which 24 of them develop lung cancer. They then want to compute the incidence rate of the process and proceed as follows: $$\text{Incidence rate} = \frac{24}{2300-24/2}$$ At first, I didn't understand why they subtracted $24/2$, but I assumed it was some correction for the fact that since those 24 persons develop the cancer during the year, their time at risk is shorter than that of the ones not developing the disease. No further information was given in the textbook itself, at least not in the problem. A quick search confirmed that I'm thinking along the correct lines. But I still don't understand the rationale for the formula. Can someone enlighten me? Also, if some references accessible to medical students could be given. I don't mind having more techical references as well.


title: Log - likelihood function, why does the summation sign vanish?

text: I have the log-likelihood function: $$l(p_i,y_i) = \sum_{i = 1}^n \left( \ln(p_i) + y_i \ln(1 - p_i) \right) $$ And I need to calculate the maximum likelihood estimator of $p_i$. When I do this, for some reason when differentiating, the summation sign vanishes. Why is this? EDIT: To calculate the maximum likelihood estimators, I would differentiate my log likelihood, equal that to 0 and solve: $$ \frac{\partial{l}}{\partial{p_i}} = \sum_{i = 1}^n \left( \frac{1}{p_i} - \frac{y_i}{1 - p_i} \right) = 0 $$ Rearranging gives me $$ \sum \frac{1}{p_i} = \sum \left( \frac{y_i}{1 - p_i} \right) $$ Oh, does this now become $$\frac{n}{p_i} = \frac{n y_i}{ 1 - pi} $$ And then I can divide through by n and rearrange to get my value for $\hat{p_i}$?


title: How can I find a measure of fit for two different sets of data?

text: I've got two different sets of data. Each set of data is a two dimensional matrix, filled with real values. What I want to do is measure how well these sets of data correlate (so that a perfect fit would have a value of zero). For example, if the two sets of data, when plotted, are each very far apart from each other - this should be not very well correlated, and have a large value associated with it. However, if the data sets are similarly shaped and positioned very similarly, the 'correlation value' should be small. The catch is that these data sets don't have matching x or y columns - in fact, each data set has a different number of elements. I was going to use the least squares method but then realized that this wouldn't work that easily in this case. Are there any well known or documented techniques for doing this? If not, any ideas how I might be able to achieve this?


title: How to edit features of a chart from SPSS in MS Word?

text: Is there any way to edit a scatterplot created in SPSS/PASW after exporting/copying and pasting it into a Word document. I am finding that when I copy and paste my graph from SPSS/PASW into Word, the font size decreases. I can enlarge the graph in Word by dragging the corners but this makes the font look odd. I want to be able to edit the font size in Word. Is there any trick that I need to learn. I am very new to SPSS/PASW.


title: How to find when a graph reaches a peak and plateaus?

text: This may sound very basic, but I have this problem: I've got a queue of data with a window size of 300. New data is added at one end, old values are removed from the other end. I expect the queue data to stay more or less consistent, e.g.: 10,12,15,10,20, then start to rise sharply: 15,10,20,22,25,26,28,30,32... all the way up to 150 or so. There the data may fluctuate a bit, then it would go down with a similar slope (120,118,116,115...) all the way down to 20 or so. I'm trying to identify the turning point in this data series programmatically, but my code detects peaks more often than I'd like to. How can I identify when the graph is rising, when it has reached a definitive turning point and when it starts to fall? Should I try to look at the rate of change of the rate of change?


title: Validate a medical test

text: I've never done something like this before, thus I do not even know where to start with calculating whatever I need. Following scenario: We measure a medical value (Glucose) with our device, and also measure a laboratory reference value to compare our value to. Now for a paper I was asked to deliver some sort of value which specifies how good our measurements are. I know of correlation as a measurement of how dependent of each other two measurements are, but I have no idea what sort of value is being asked for in this sort of tests? I read something about Level of Confidence, which sounds good, but I don't understand how I would apply it here. The data I have is (basically) two arrays, each with 30 floating values. What do I do from here?


title: Prediction intervals for a Poisson count (in R)

text: Is there a convenient way in R to calculate a probability prediction interval for a count sampled from Poisson($\lambda$), for known lambda? That is, for a given $\alpha$, find $a$ and $b$ such that $P(a\le X\le b)\ge1-\alpha$, where $X\sim\text{Pois}(\lambda)$. Is there a generic way to do this in R for any of its built-in distributions?


title: Asymmetric S-shaped function mapping interval $[0, 1]$ to interval $[0, 1]$

text: In the literature, is there any asymmetric $S$-shaped function that maps the interval $[0, 1]$ to interval $[0, 1]$? Unfortunately I can't post figure so I just describe what I mean in text. The function I want should be monotonic increasing. In the meanwhile, for a real number $c < 0.5$, the function $$0 \mapsto 0, \quad c \mapsto 0.5, \quad 1 \mapsto 1.$$ Moreover, the function is smooth on its domain, concave on interval $[0, c]$ and convex on interval $[c, 1]$. By the way, can it be seen as a sort of link-function?


title: Variable selection for linear regression using robust or least squares estimation

text: I have a data set consisting of one continuous response variable and about 70 predictors. Using this data, I want to construct a linear regression model. However, I don't know what predictors are worth including in the model, so I'll need to utilize a variable selection method that will allow me to isolate specific response variables. Unfortunately, I've noticed that my data violates a number of assumptions associated with linear regression. Therefore, I'll need to utilize a different estimation than OLS, such as robust or least squares estimation. When running a linear regression with a different estimation method than least squares, how does one utilize a variable selection method such as stepwise? How can this be implemented in R?


title: Why is there an R^2 value (and what is determining it) when lm has no variance in the predicted value?

text: Consider the following R code:               example <- function(n) {         X <- 1:n         Y <- rep(1,n)         return(lm(Y~X))     }     #(2.13.0, i386-pc-mingw32)     summary(example(7))    #R^2 = .1963     summary(example(62))   #R^2 = .4529     summary(example(4540)) #R^2 = .7832     summary(example(104))) #R^2 = 0     #I did a search for n 6:10000, the result for R^2 is NaN for     #n = 2, 4, 16, 64, 256, 1024, 2085 (not a typo), 4096, 6175 (not a typo), and 8340 (not a typo)      Looking at http://svn.r-project.org/R/trunk/src/appl/dqrls.f) did not help me understand what is going on, because I do not know Fortran. In another question it was answered that floating point machine tolerance errors were are to blame for coefficients for X that are close to, but not quite 0. $R^2$ is greater when the value for `coef(example(n))["X"]` is closer to 0. But...   1. Why is there an $R^2$ value at all?    2. What (specifically) is determining it?   3. Why the seemingly orderly progression of `NaN` results?   4. Why the violations of that progression?   5. What of this is 'expected' behavior?


title: Significance of difference using a distance matrix

text: So my data consists of a distance matrix for some points, and a table classifying the points as either red or green. I want to know if there is any difference between the red and green points. I did UPGMA clustering and they don't look different (doesn't seem to be any tendency for red or green to cluster together) but I don't know how to express the significance of this. I guess the question would be something like "what is the chance of getting a distance matrix like this, if the points were all randomly distributed".


title: How to analyse seasonal variation in stroke onset?

text: I am dealing with the topic of seasonal variation and stroke. I would like to do time series analysis.For time variable I have stroke onset month and all the other variables are categorical in nature. How should I analyze these data?


title: Difference between ep-SVR and nu-SVR (and least squares SVR)

text: I am trying to find out which SVR is suited for what kind of data. I found out 4 types of SVRs: epsilon, nu, least squares and linear. I understand linear SVR is more or less like lasso with L1 Reg. But what is the difference between the remaining 3 techniques?


title: Determining the uncertainty of an exponential fit

text: My problem should probably be built up from the beginning, so lets start there. I performed a certain experiment 25 times. Every time, the experiment consists of 5000 measurements, and each measurement returns either 1 ('yes') or 0 ('no'). I know that each time I measure 1, the probability that this measurement is correct is $F_1$, while if I measure 0 I know that it's correct with probability $F_0$. Now, subsequently I want to compute the autocorrelation function of these measurements I use the following formula ![enter image description here](http://i.stack.imgur.com/G3L64.png) with a slightly different normalization, but in principle that should not matter, because what I want to do is fit this autocorrelation to an exponential decay and find the decay time. To do so, I thought I had two options: Either I fit the 25 datasets separately, and then find the mean of the decay time $t_1$ and maybe say something about its variance, or I can find the mean autocorrelation of the 25 datasets, and fit that one. I've decided that that is probably the best option, as the individual datasets can be quite noisy with strange fluctuations, while the mean looks much cleaner. But now I'm stuck wondering what I should be doing to find the uncertainty in the fit. Should I have introduced some sort of uncertainty due to the imperfect measurements ($F_0$ and $F_1$), or should I have introduced some sort of standard error of the mean for the 'mean autocorrelation' data? That last one seems plausible, as of course there is also a standard deviation in the 25 datasets, but then I get a little confused as to how I should do that. The formula for the standard error of the mean seems to be $\frac{\sigma}{\sqrt{n}}$, but this is a bit unfair as in my 5000 measurements, the autocorrelation drops to 0 after around 50 measurements. For the fit I therefore also only use the first ~100 points of the autocorrelation, as there's no real need to fit the next 4000 points that are all pretty much equal to 0. So should I instead only calculate the mean and standard error of the first 100 points? Would that be a fair way of finding the uncertainty? The subsequent fitting method will already give error bars for the $t_1$, but this of course heavily depends on the uncertainty in the data. A final sort of sidetrail, what kind of quantitative method of finding how good I fit should I use for exponential decay? Reduced Chi square, Adjusted R squared, or something else entirely?


title: What's the problem with model identifiability?

text: I understand that in a decision perspective, identifiability of a model is needed to ensure the convergence (with increasing number of observations) of the parameters to estimate through a single value. But, if the non- identifiability of a given model is not a modeling artifacts but clearly characterises some "inaccessible knowledge" about the system under study, is it valid to perform bayesian inference on a non-identifiable model ? Here is a simple example. $$ x_i =t a y_i + \epsilon_i $$ with $(\epsilon_i)$ iid $$ \epsilon_i \sim N(0,1) $$ and an informative prior for $t$: $$ t\sim N(1,0.1) $$ and a non-informative prior for $a$ (let says, that one chooses a uniform...) $$ a \sim U(0,1000) $$ One observes $(x_i)$ and $(y_i)$ are exogenous parameters and one wants to compute : $$ p(a | (x_i); (y_i)) $$ As I understand it, the model is not identifiable as all the densities $p((x_i) | a,t;(y_i))$ described by the pairs $(a,t)$ such that $a.t=k$ ($k \in R$) are the same. Obviously in such a case the choice of $p(t)$ has a strong implication but if it is physically supported, I see no reason to invalidate the meaning of an HPD interval obtained from such a non-identifiable model. On the other hand, I do not manage to find any reference about that... so thanks for your expertise.


title: Using weights on dummy variables in R

text: I am doing a regression analysis using dummy variables, but I need to include weights. I am using data from an opinion poll, and I am trying to make sure that the sample is regionally representative. So I would like to apply weights for responses from certain regions. How do I do this with dummy variables? Can I simply apply to weight to the dummy before putting it into the regression? Any guidance would be appreciated. thanks


title: Formation of the test statistic in one-way ANOVA

text: From nist for the test statistic in one-way ANOVA: > The mean squares are formed by dividing the sum of squares by the associated > degrees of freedom. > > Let $N = \sum_i n_i$. Then, the degrees of freedom for treatment, $DFT = k - > 1$, and the degrees of freedom for error, $DFE = N - k$. > > The corresponding mean squares are: > > $MST = SST / DFT$ > > $MSE = SSE / DFE$ I was wondering why the degrees of freedom for the two sums are given as $k-1$ and $N-k$, where $N$ is the sample size and $k$ is the number of groups? For the definitions of the two sums $SST$ (sum of squares of treatments) and $SSE$ (sum of squares of errors), see another page. So dividing a sum of squares by its degree of freedom isn't "mean square", isn't it? Why does it say "the mean squares are formed by dividing the sum of squares by the associated degrees of freedom"? Thanks and regards!


title: Feature Selection for look alike modeling using k-NN

text: I have a list of items and various parameters for each items. For each item on my list i need to identify items which are similar to the item from my whole population . I am planning on using K-NN algorithm for the same. However, however i am not sure how to pick up the relevant parameters. Can anyone guide me on what would be an ideal technique for feature selection?


title: How to compare two different clustering approaches?

text: I'm making a project connected with identifying the dynamics of sales. My database concerns 26 weeks (so equally in 26 time-series observations) after launching the product, 126 time-series=126 different products. I used clustering to group products with similar dynamics of sales. I applied two different approaches to pam() = k-medoids clustering. First one groupes time-series - 3 clusters, second one groupes S-curve parameters, counted before by curve fitting to each time-series - 6 clusters. Now I would like to compare these two methods, but have no idea how to compare the quality of two clustering methods if I have different number of clusters. I would be grateful for any suggestions.


title: Confidence intervals for glmer() from lme4

text: I'm running a mixed model on some data. I want to calculate confidence intervals for my model. For this I have adapted the following code section from Predictions and/or confidence (or prediction) intervals on predictions (lme4). The problem is that I have a hard time understanding what the code actually does, but so far it's the only way I have found to calculate some reasonable confidence intervals for my result. So my question is: what are `pvar1` and `tvar1`? I'm guessing it is something like $\text{SE}^2$, but I can't quite see it, and I can't describe their difference, or their normal statistical term. I'm guessing that the difference is one is CI while the other is PI? Further when the code then calculates the confidence intervals it says `prediction +/- 2*sqrt(pvar1)`. Shouldn't it be the t.crit value or 1.96 for a normal distribution?               library(lme4)     library(ggplot2) # Plotting     data("Orthodont",package="MEMSS")     fm1 <- lmer(         formula = distance ~ age*Sex + (age|Subject)         , data = Orthodont     )     newdat <- expand.grid(         age=c(8,10,12,14)         , Sex=c("Male","Female")         , distance = 0     )     mm <- model.matrix(terms(fm1),newdat)     newdat$distance <- mm %*% fixef(fm1)     pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))     tvar1 <- pvar1+VarCorr(fm1)$Subject[1]  ## must be adapted for more complex models     newdat <- data.frame(         newdat         , plo = newdat$distance-2*sqrt(pvar1)         , phi = newdat$distance+2*sqrt(pvar1)         , tlo = newdat$distance-2*sqrt(tvar1)         , thi = newdat$distance+2*sqrt(tvar1)     )     #plot confidence     g0 <- ggplot(newdat, aes(x=age, y=distance, colour=Sex))+geom_point()     g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+         opts(title="CI based on fixed-effects uncertainty ONLY")     #plot prediction     g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+         opts(title="CI based on FE uncertainty + RE variance")      Thank you for your help.


title: Difference test between 2 groups where both groups contain depended and independent observations

text: Is there a statistical test to check the difference between groups (for example a paired sample proportion test), but where some of the observations of the two groups are dependent and some of the observations are independent? **For example:** **Context** : The observations of the groups are named with a letter and can have a value of 1 and 0 ![enter image description here](http://i.stack.imgur.com/taXkV.png) **dependent observations:** A, D, G & Q **Independent observations** : B, C, F, E, M, N , O , P , R & S **Hypothesis I want to test:** There’s a significant difference of the proportion of 1’s between the 2 groups. **Problem:** I can’t use the dependent proportion test because some of the observations in the groups are independent. And I can’t use a independent t-test because some of the observations are independent. Thanks in advance.


title: Checking if two Poisson samples have the same mean

text: This is an elementary question, but I wasn't able to find the answer. I have two measurements: n1 events in time t1 and n2 events in time t2, both produced (say) by Poisson processes with possibly-different lambda values. This is actually from a news article, which essentially claims that since $n_1/t_1\neq n_2/t_2$ that the two are different, but I'm not sure that the claim is valid. Suppose that the time periods were not chosen maliciously (to maximize the events in one or the other). Can I just do a _t_ -test, or would that not be appropriate? The number of events is too small for me to comfortably call the distributions approximately normal.


title: Yates' correction for continuity only for 2X2?

text: Is the Yates' correction for continuity used only for 2X2 matrices?


title: Validate web a/b tests by re-running an experiment - is this valid?

text: A webinar the other day by an a/b testing company had their resident "Data Scientist" explain that you should validate your results by re-running the experiment. The premise was, if you select 95% confidence, there is 5% (1/20) chance of a false positive. If you re-run your experiment with the same constraints, now there is a 1/400 (I'm assuming they determined this as 0.05^2 = 1/400) Is this a valid statement? (ie, "run twice, two statistical significance wins = 1/400 probability of false positive")? Would it have been a better approach to increase your significance level? From a business standpoint, the concern I have is by re-running the experiment, you are exposing more users to an inferior page (treatment), and thus losing out on potential sales.


title: Determining if change in a time series is statistically significant

text: I have the total number of calls received each week and have plotted them on a chart, going back nearly 3 years. By eye it seems that there was a massive drop over Christmas, that doesn't seem to have recovered, it seems that there has been a step change in requests. Is there a test I can do that can quantify this difference? Cheers Ben


title: How to handle variables with low correlation but high loadings in factor analysis

text: I am doing factor analysis to check the factorial validity of a 14-items scale with four subscales. Two items have low (less than 0.3) correlations with other items in the subscale to which they belong. However, they load considerably (greater than 0.32) on the factor (latent variable) to which they correspond. Under such condition, should the two variables be dropped from the scale on grounds of low correlations or be retained within the scale on grounds of the acceptable loadings? Update: After going further through the analysis--removing one item on substantive ground and another item due to low loading (in the range of .2)--I found the correlation matrix depicted below. ![enter image description here](http://i.stack.imgur.com/HFPT5.png) The highlighted correlations indicate the correlations of items with other items in the same subscale. However, some items also correlate with items in other subscales considerably (>0.3). Yet, the four factors themselves are inter-correlated which I think is the reason why items also show non- negligible correlations with other items of a different factor. Below is the inter-correlation of factors. ![enter image description here](http://i.stack.imgur.com/A6S76.png) However, as you see hereunder, the loading of each item on the factor to which it belongs is substantial and no cross-loading. ![enter image description here](http://i.stack.imgur.com/1TVxg.png) Given these results, how should I go about interpreting my subscales? May I be justified to consider the sub-scales as being suited to measure the underlying constructs/factors of interest?


title: What is the mixture of distribution hypothesis?

text: I keep seeing this term crop up in financial math papers but am unable to find a definition/explanation without purchasing whole articles! Can someone shed some light?


title: PCA figure formatting options in R

text: I've completed PCA analysis, in R with VEGAN package, of some ecological data on tree health. There are 80 trees total (so, 80 'sites') divided into four treatment categories. I've got the data plotted with color coded points-- colors according to the treatment groups. Rather than plotting individual sites/trees on PCA biplot, I'd like to make something like a box-and-whisker plot that has four 'crosses' that show the centroid for each group and the SE in both PCA dimensions. I've seen figures like this in papers, but I can't seem to find an R script for plotting this way. Any suggestions? I guess an alternative would be to just take the site scores and manually find the means and SE's and create my own plot, but I'd rather find a script for it, if possible.


title: Variational inference engines

text: After doing some research on the topic, I have noticed a surprising deficit of inference packages and libraries that rely on message-passing or optimization methods for Python and R. To the best of my knowledge, these methods are extremely useful. For example, for a Bayes Network (directed, acyclic) belief-propagation alone should be able to give exact answers. However, most inference software that is available online (e.g. STAN, BUGS, PyMC) rely on MCMC methods. In the Python case, to the best of my knowledge, neither PyMC, scikit-learn or statsmodels include variational inference algorithms such as belief propagation, message-passing methods or any of their variants. Why is that? Are these methods less used in practice because they are seen not as powerful or generic as their MCMC counterparts? or Is it simply a matter of lack of manpower and time?


title: What is a mixing process?

text: What does this mean? > Asset prices follow a mixture of normal distributions with a mixing process > dependent on the unobservable information arrival process.


title: CIs for proportions in Stata: ci vs. proportion command

text: I'm performing some basic analysis on a large data set (n-size after restrictions about 4,000) to supplement some qualitative historical research. The goal is to produce a variety of confidence intervals for proportions, and I've run the analysis in Stata using both the `ci` command (with `aweights`) and the `proportion` command (with `pweights`). The means generated by `ci` are more or less identical to the proportions generated by `proportion`, but there is more noticeable (though not very substantial) variation in the confidence intervals. I'm wondering if there is a strong theoretical reason to opt for either the `ci` or the `proportion` command here? I'm guessing the latter given the sample size.


title: Why do I get more than 25 symbols in plot in 25-fold cross-validation

text: I have a group of observations, 25, And I want to perform a 25-fold cross- validation (leave-one-out). Here is my code,               library(MASS)     min.m <- lm(formula= Mass ~ 1, data=mean.mydata)     fsel.m <- step(min.m, direction = "both", scope=sel.fmla, steps = 5)          fsel.fmla <- as.formula(object=fsel.m)          # Leave-one-out rule validation     library(DAAG)     cross.m <- suppressWarnings(cv.lm(df=mean.mydata, fsel.fmla, m=25, legend.pos=FALSE))      This is my plot, ![enter image description here](http://i.stack.imgur.com/mAu1J.png) Can anyone explain me why there is more than 25 symbols? I noticed there is 2 of each ones, does it correspond to cross-validation results (cvpred) and predicted data?


title: Why are the geometric distribution and hypergeometric distribution called as such?

text: Why are the geometric distribution and hypergeometric distribution called "geometric" and "hypergoemetric" respectively? Is it because their pmfs take some special form? Thanks!


title: Binomial-like distribution

text: Let's say that $F(k; n, p) = \Pr(k \leq X)$ is the cumulative binomial distribution for $X \sim B(n, p)$. I am trying to solve for the smallest $i$ such that $F(i; n + i, p) > \alpha$ for some $\alpha$. I would use the inverse binomial distribution function to find this, except the $n$ parameter changes. Does anyone know if this follows any known distribution?


title: Is there a measure of association for a nominal DV and an interval IV?

text: As part of a larger analysis, I'm interested in examining the relationship between an interval level variable and a four category nominal variable. In this case, the interval variable is the independent variable and the four category nominal variable is the dependent variable. I'm interested in whether there is a statistically significant relationship between the two variables, as well as the strength of the relationship if one exists. Does anyone know whether a measure of association exists for this situation? I'm aware of eta, but my understanding is that eta is appropriate when the IV is nominal and the DV is interval. I'm also aware of the point bi-serial correlation, but my understanding is that the nominal variable must be binary. Any help is greatly appreciated.


title: GARCH(1,1) Implementation question

text: Could someone shed some light into the implementation of the GARCH(1,1) model contained in page6 of the following document? http://www- stat.wharton.upenn.edu/~steele/Courses/956/RResources/GarchAndR/WurtzEtAlGarch.pdf. In particular, I am looking at the creating the likelihood function in R. To begin with, the GARCH(1,1) model is define to be $\sigma _t ^2 = \omega + \alpha z^2 _{t-1} + \beta \sigma ^2_{t-1}$. By repeated substitution, the model can be written as $\sigma _t ^2 = \omega + \alpha z^2 _{t-1} + \beta(\omega + \alpha z^2 _{t-2}) + \ldots \beta ^{t-2}(\omega + \alpha z^2 _1) + \beta ^{t-1}\sigma_1 ^2$. Now, based from my understanding of the use of the function filter in R, the line of code `h = filter(e, beta, "r", init = Mean)`, seems to suggest that $\sigma ^2_t \equiv \omega + \alpha z^2 _{t-1} + \beta(\omega + \alpha z^2 _{t-2}) + \ldots \beta ^{t-2}(\omega + \alpha z^2 _1) + \beta ^{t-1}(\omega + \alpha \mathbb{E}[z^2]) + \underbrace{\mathbb{E}[z^2]}_{=:\text{init}}$, which is somewhat inconsistent with the previous equation. I am not aware of any theory which suggest this relationship or at least I am not sure why the code initialises the historical variance of the data outside the recursive equation. I hope somebody can help me address this issue for me if I have misunderstood it.


title: Is it possible to determine the probability of NOT reaching a state in a FSA ?

text: Suppose I have a Finite State Automata with various states and probabilities of state transitions. Does mathematics exist to determine the probability of NOT reaching a state in the FSA given some n number of state changes?


title: computing mean and sd of d-prime from mean and sd of hits and false alarms

text: I'm working on a meta-analysis. Some studies report d', others report accuracy (hits and false alarms). I would like to transform all the data into the same measure, d' (also c). However, for many studies, I only have the mean and sd of the accuracy measures. If I compute d' from the mean hit and false alarms, the outcome is not the same as the mean d', but it is at least closely correlated with mean d'. However, I cannot figure out how to figure out the variance of the d' based on the mean and sd of the accuracy data. Is this possible, or am I going to have to request individual subject data?


title: Create age and sexmatched pairs to balance Cox regression further (updated)

text: I analyze ethnic differences in risk of cardiovascular events (CVD) in a cohort study of patients with coronary heart disease. It is known that immigrants have higher risk of CVD and I intend to show this via Cox regression. However, my results fail to prove what is both plausible and well known, i.e immigrants do not have significantly higher risk. I suspect the reason for this is that groups are not comparable, not even with Cox model, because immigrants are 10–20 years younger at baseline (age is an important risk factor). There might be insufficient number of individuals with overlapping age to compare immigrants to natives. So I thought the comparison would be more appropriate if each immigrant would be matched to a sex and age matched control. I used the MatchIt package (I have some basic understanding of prop scores and the fact that they are primarily intended for treatment exposures). I reasoned that, as treatment is an exposure, so is immigrant status, which is dichotomized (yes/no immigrant status). I obtained (adjustments for age, sex and duration of coronary heart disease) the probability of being native. I used the weights in coxph function as follows:               coxph(Surv(start, stop, event==1) ~  var1 + var2 + var3 + var4, data=dataset, weights=weights)      Results are now more plausible. Immigrants are now at higher risk. Was this correct or just a incorrect effort to prove an hypothesis? Adam


title: tease out significance

text: Say I have search data like this `AvgCost QualityScore SearchShare` `3.12 6 0.6364` Where AvgCost is a continuous numerical variable, qualityScore is a categorical variable with values 1-10, and SearchShare is a percentage...I am wondering how to tease out the effect Avgcost and QualityScore have as well as whether they are too related to be both included. I'm also wondering what kinds of regressions besides simple linear regressions can be run on data sets where I have these sorts of variables i.e. categorical and numerical explanatory variables and a percentage response variable. I have done some obvious things, like check `cor()` between qualityScore and AvgCost as well as do a simple linear regression. I looked at the `relaimpo` package, not sure if that's what I need. Anyway, thanks very much for any help!


title: Sample Range to Estimate StDev

text: Accidentally, I found (by simulation) that the range of a normally distributed sample is related to the estimator of the standard deviation: $$X \sim \mathcal{N}\left(\mu, \sigma\right)$$ $$s(X) = \frac{1}{3}\left(\max(X) - \min(X)\right)$$ Why is that? I can't find any references. Is it related to the three sigma rule?


title: The first principal component becomes irrelevant

text: I did run PCA on 17 quantitative variables in order to obtain a smaller set of variables that is principal components to be used in supervised machine learning for classifying instances into two classes. After PCA the PC1 accounts for 31% of the variance in the data, PC2 accounts for 17%, PC3 accounts for 10%, PC4 accounts for 8%, PC5 accounts for 7% and PC6 accounts for 6%. However, when I look at mean differences between PCs among the two classes, surprisingly, PC1 is not a good discriminator between the two classes. Remaining PCs are good discriminators. In addition, PC1 becomes irrelevant when used in a decision tree which means that after tree pruning it is not even present in the tree. The tree consists of PC2-PC6. Is there any explanation for this phenomenon? Can it be something wrong with the derived variables? Thanks!


title: How to mine the "customers who bought X also bought Y" type of correlation?

text: I am looking for a statistical model to detect grouping patterns among a population X of n element with respect to their association with k elements of population Y. In modern marketing, an example of that would be "customers who bought X also bought Y", such as on Amazon etc.               X (person)  Y (purchase)     ------------------------     Jimmy       Lure     Joe         Fishing Rod     Jimmy       Fishing Rod     Sue         Sewing Machine     Sue         Thread     Patty       Sewing Machine     Jack        Tennis Balls      Above is just a very simple example. But imagine that 80% of people who bought a `Fishing Rod` also bought `Lure` and similar for `Sewing Machine` vs `Thread`. While the pattern may be obvious to the naked eye, I am looking for a more deterministic algorithm to detect this seemingly simple pattern. Imagine that this is a simple two dimensional array. Could you produce a pseudo code (or Java, Python, C or whatever suits you) algorithm that I could use to determine the most closely related elements in X based on their association with elements from Y?


title: Linearized Ridge Regression Estimator Under the Mean Squared Error Criterion in a Linear Regression Model

text: I have been implementing this paper: http://www.tandfonline.com/doi/abs/10.1080/03610918.2011.575506#.UrijOLS_UhI I can send a copy should anybody require it. I think that my implementation is correct but I get significantly different answers to those given in the paper, when I apply the LRRE to the Hald data set. (Note that in the paper the OLS solutions for the Hald data set are given and these match exactly so I believe the data is correct.) My results   B0 0.0470   B1 2.1635   B2 1.1560   B3 0.7451   B4 0.4931 Results in paper   B0 0.1228   B1 2.1630   B2 1.1552   B3 0.7444   B4 0.4922 Has anyone else implemented this paper? Did you also get different results to those quoted by the author? I'm happy to share my code. The data is as follows Here is the Hald data set:                78.5   1   7   26  6   60      74.3   1   1   29  15  52      104.3  1   11  56  8   20      87.6   1   11  31  8   47      95.9   1   7   52  6   33      109.2  1   11  55  9   22      102.7  1   3   71  17  6      72.5   1   1   31  22  44      93.1   1   2   54  18  22      115.9  1   21  47  4   26      83.8   1   1   40  23  34      113.3  1   11  66  9   12      109.4  1   10  68  8   12


title: Test whether a set of samples follow laplace distribution

text: I have a set of samples from a random variable, is there any method so that I can confirm/check to know whether they follow a laplace distribution in an statistically rigorous way? I am very new to statistics and would really be grateful if you could at least point me towards some litrature. Thanks


title: Bootstrapping Average Treatment Effects and Confidence Intervals for multinomial logit models

text: I have a multinomial logit model with three outcome categories and a variety of (primarily binary) independent variables. How would I go about bootstrapping the average treatment effect and its confidence interval of a binary variable? Right now my procedure is as follows: After running the full model as a multinomial logit I   1. Sample with replacement to create temporary data of same length as original dataset.   2. Reestimate the model using temporary data   3. set set indep. variable to 0 in temporary data and predict probabilities of outcomes (using reestimated model)   4. set set indep. variable to 1 in temporary data and predict probabilities of outcomes (using reestimated model)   5. Take the mean of the average difference in predicted probabilities   6. Save this mean difference in pred probabilities   7. Repeat this process a lot of times   8. Take mean of all the mean differences pred. probabilities for average treatment effect, Take 0.25 and 0.975 quantile for a 95% confidence interval around the mean.  The issue I am having right now is, that this seems to create very large confidence intervals. This has created a situation where a lot of the variables that are significant in the regression output no longer have an effect that is differentiable from 0. Two questions:   1. Is this bootstrap procedure sensible?   2. If so, what explains the difference between the regression outputs and the bootstrap results?


title: Sentiment Analysis with respect to subject

text: I'm familiar with the bag of words/Naive bayes sentiment analysis for text (e.g. http://streamhacker.com/2010/05/10/text-classification-sentiment- analysis-naive-bayes-classifier/) I was curious to see if anyone here knows anything about computing sentiment of a sentence with respect to a subject. For example, with respect to the word 'Joe', the sentence "I am happy" is neutral sentiment whereas "Joe is happy" is positive sentiment Do the traditional bag of words features work for this type of sentiment classification? or is a more sophisticated feature set necessary? Thanks!


title: Can we never have a normal distribution?

text: Since the population variance can never be estimated, will we always have a t-distribution instead of normal?