title: Efficient generation of chi random variables

text: I need to generate random variables generated from a chi distribution ( **not** chi-squared!). There doesn't seem to be standard mechanism in C++ in (for example) Boost::Random and hence I am looking for an alternative implementation implemented in C++ (or an easily compatible language such as C or FORTRAN). There is one paper in ACM TOMS on the topic. Or, should I simply take the square root of a chi-squared random variable, of which there are several implementations? Which is likely to be fastest - unfortunately, I need to draw a lot of them. * * * **Update:** In response to the below comments, these are some notes for my specific application:   * I am looking for a C/C++ compatible implementation, or sufficient detail to program one myself.   * The application is CPU sensitive, and is run on a modest embedded platform (400MHz with a scalar double FPU and limited cache, but no sqrt(), sin(), exp() or similar instructions). Since actual run-time is not easy to predict _a priori_ , I am looking for algorithms that are _amongst_ the fastest and known to perform well. I can evaluate each of these on my own hardware.   * I already have facilities for generating uniform and normal random numbers. I can get chi-squared if necessary.   * DoF is always an integer, usually between 15 and 39. The answer I accept will have the following content:   * Which methods are considered "fast", and how are they implemented?   * Do they otherwise have known disadvantages (e.g. numerical stability)?   * Are there any methods that should be specifically avoided?


title: Generating and verifying uniformly distrubuted orthogonal matrices

text: I would like to generate a uniformly distributed $n \times n$ orthogonal matrix. There seems to be several such methods; see this question and the oft- cited paper by Stewart. Using a QR decompostion of a matrix of normally- distributed variables seems to be another popular method. However, I am being careful, because apparently it is not hard to generate a subtly-wrong distribution. This paper by Mezzadri demonstrates how it can be wrong for a random _unitary_ matrix. So, I would like to ask:   * What implementation issues do I need to be aware of when applying either the QR or the Stewart method?   * Are these indeed the fastest methods? (I want to generate $15 \times 15$ matrices on an embedded platform, so both CPU and cache are of the essence)   * How to I verify that the generated distribution is correct? The last question is important - How do I verify that distribution is correct? What do I need to look for? How can I implement such a test, even if visually for (say) a $2 \times 2$ or $3 \times 3$ matrix?


title: Comparing SEM software apps: AMOS vs Stata vs Mplus vs R ---- opinions please?

text: I have AMOS and Stata, and I have a colleague who has Mplus (and I'd rather not ask the colleague to run my data for me). I recently discovered that Stata has SEM capabilities. I looked online, but I can't find any comparison of Stata vs AMOS vs Mplus... What are your thoughts on comparing the 3? I'm particularly interested in opinions about SEM in Stata because I can't find anything online discussing how good it is.


title: What is maximum number of variables that we can use to run a cluster analysis?

text: To do one clustering analysis, the model that I developed contains some 30 variables. I need to run this clustering for some 2-3 million data points. I need to know whether number of variables that I intend to use (30) is just too many? Would my result be impacted 'coz of large number of variables? Also there would be issues with processing limitations of tool (spotfire) that I am using for this exercise. But that is secondary concern at this moment.


title: How to find quantiles for multivariate data using R?

text: Quantile for a single variable is easy to implement in R. However, it is not an easy task to quantile for multivariate data. There are several papers have been proposed to quantile for multivariate data such as Chaudhuri, P (1996), On a geometric notion of quantiles for multivariate data, _Journal of the American Statistical Association_ 91, 862-872. The existing approaching is not easy to implement. Is there any R package can do quantiles for multivariate data? many thanks in advance


title: Bootstrap distribution of a breakpoint not centered around its parameter estimate

text: I estimated the breakpoint for a piecewise linear model (one intercept, two slopes that meet at the breakpoint) by minimizing the deviance using an optimizer (`optim` in R). I noticed that the estimate was near 2.0, 4.0, and 6.0, depending on the initial conditions, indicating there were local optima. So I used simulated annealing to estimate the breakpoint (`optim(..., method="SANN", ...)`), and became reasonably confidence that I had found the globally best estimate at bp = 6.0. I just finished bootstrapping the breakpoint estimates (my data are nested within subjects, and I had to apply simulated annealing to each bootstrap estimate, so the whole thing took a couple days), and found that: 1) The bootstrap distribution is bimodal, with modes at 2.0 and 6.0. 2) The mean of the bootstrap distribution is 4.0, and the median is 3.5, whereas the parameter estimate was 6.0. Can someone help me to understand my results? I want to say that the original estimate of 6.0 is an artifact of my sample, and that the bootstrap mean of 4.0 is a better estimate of the breakpoint, but I'm not sure if this is the right interpretation. Also, placing the breakpoint at 4.0 produces conditional means that agree almost exactly with a LOESS curve using default values (a breakpoint at 6.0 does not agree as well). \-- Edit: Adding a graph. ![Red line = original fit with bp = 6.0. Black solid = fit with bp = 4.6. Black dotted = LOESS. Open circles = means. Models were fit to full data set \(too many data points to show, and yes, I tried high density plotting methods as well\), not means.](http://i.stack.imgur.com/vnkLS.jpg) Figure notes: Red line = original fit with bp = 6.0. Black solid = fit with bp = 4.6. Black dotted = LOESS. Open circles = means. Models were fit to full data set (too many data points to show, and yes, I tried high density plotting methods as well), not means.


title: Low intensity Poisson estimation

text: I have a collection of Poisson processes each with an unknown $\lambda$. I would like to estimate $\lambda$ for each process. for each process I could take either the total number of event over the total time or the inverse of the mean of the waiting times. Given that the dataset is quite small, with many of the processes having less than 5 total events, is is better to use the estimate based on the waiting times or the total count? The replies to this question suggest that the count over the time is the best estimator, but no mention of the waiting times.


title: How to resample matrices to test for the robustness of their correlation?

text: I have several populations where I have morphology and diet for each individual. I am interested in the correlation between diet and morphological distances. However the number of individuals in each population ranges from 22 to 80 individuals. I have looked at the correlation diet-morphology for each population and (not surprisingly) the correlation coefficientis is highly correlated with the number of individuals per population. I would like to resample (without replacement) the populations with 60-80 individuals and get random samplings of 30 individuals (1000 times). I would like to get a correlation coefficient distribution against which to test the original value of the correlation. I guess it is possible do this in R, however I have never written a script in R and I am not familiar with resampling techniques at all. Any help with coding will be greatly appreciated Thank you Camille


title: Find the effect of a attribute value on an outcome by eliminating confounding values

text: I have a series of lets say five attributes. The first attribute is called diagnosis code 1, the second diagnosis code 2 etc. The values are codes which represent diseases. In other words, each observation, can have up to five DIFFERENT diagnoses. I want to know the impact of one given diagnosis code on a specific outcome (Length of Stay for example), but in order to do this, I have to eliminate the confounding effect of the coexisting diseases to the length of stay. I am sort of stuck on how I can do this. Almost all patients (observations) have data in all their five diagnosis slots. This is going to be part of a software, but I would first like to understand the overall approach of how to solve this, and then we will think about the programming part. thanks!


title: One sample Hoteling's T2

text: I used one sample Hotelling's T2 statistic to check if the nutrient intake of 75 men meet the recommended daily nutrient intake value. I got an F value close to 58000. Is such a large value permissible? I checked the multivariate normality of men using MVN package of R software. Now, is there anything wrong if I compare the nutrient intake of men with that of women using Mann Whitney U test as the nutrient intake of women was not found to be normal.


title: Poisson distribution confidence interval

text: Edit : this is the data I used for the first part of the problem : \begin{matrix} Rocks & 0& 1& 2& 3& 4& 5& 6 & 7\\\ samples & 12& 27& 28& 19& 8& 3& 1 & 1\\\ \end{matrix} I did hypothesis to test if it follows poisson or binomial and concluded it could follow both. ## Part 2 Let p be the probability that a sample taken on a land contains strictly more than one rock. A) What is the p value if we suppose the following : finding rock in samples follows a Poisson P(2) with $\lambda = 2$ ? For this part I just do Poisson average is 2 and random is 1 : $0.594$ * * * B) What is the _empirical proportion_ $\hat p$ of the samples that contains $ > 1$ rock ? For this part, I have, using the continuity correction ( +0.5) : $P(X > 1 - 0.5 - 2) / \sqrt{2} = p(z > 0.5 - 2) / 1.41 = -1.0638$ $ = p(z < 1.0638)$ my answer : everything below 1.06 standard deviation left from the mean. Am I right with my interpretation of _empirical proportion ?_ * * * C) Give a 98 % confidence interval for $p$ I found an interval using this calculator but I'm not sure if it's right to use 1 as my observed events value. $[0.01, 6.64]$ also found this is the formula to use $\lambda \pm 1.96\sqrt{\lambda / n}$ but in my case it's going to be $\lambda \pm 1.96 * \sqrt{1}$ ? * * * D) **Test at $\alpha = 5%$ the hypothesis $p = 0.5$ versus $ p > 0.5$** I don't know where to start this part.


title: Is there a significantly different number of people between groups?

text: How do I estimate whether a given number of people across groups (say, 72 people in group 1, 95 in group 2 and 70 in group 3) significantly differs from a random distribution of people across groups? That is, which statistical test do I have to apply?


title: K-means cluster analysis with K=2 as a binary classifier

text: I used two variables, height and weight, and using K-means cluster analysis with $K=2$, two clusters were obtained. I used $K=2$, as the observations either belong to men or women. I then compared the obtained clusters with the real classification. I observed that K-means did pretty well. Does this sound logical?


title: How to calculate prediction intervals in major axis regression using R

text: Need help in calculation of prediction intervals in major axis regression (MA). I'm using 'lmodel2' package for calculation of the MA, but I don't understand how to calculate prediction intervals using this package. For example I used a data that comes with the package: > library(lmodel2) > > data(mod2ex1) > > Ex1.res <\- lmodel2(Predicted_by_model ~ Survival, data=mod2ex1, nperm=99) > > Ex1.res Need help in calculation of prediction intervals in major axis regression (MA). I'm using 'lmodel2' package for calculation of the MA, but I don't understand how to calculate prediction intervals using this package. for example I used a data that comes with the package: `>library(lmodel2)` `>data(mod2ex1)` `>Ex1.res <- lmodel2(Predicted_by_model ~ Survival, data=mod2ex1, nperm=99)` `>Ex1.res` Here I stack. Actually I need something like a `predict()` function that calculate prediction intervals for linear regression model `lm()`. For example: `>lm_data = data.frame( Predicted_by_model=mod2ex1$Predicted_by_model,Survival=mod2ex1$Survival)` `> mod1 <- lm(lm_data)` `> #I want to know the prediction interval of variable 'Predicted_by_model' when Survival=2` `> newdata = data.frame(Survival=2)` `> predict(mod1, newdata2, interval="predict")` > >       fit     lwr      upr >   > > 1 2.000688 1.73913 2.262245 Thanks, Ilia


title: How to fit mixture model for clustering

text: I have two variables - X and Y and I need to make cluster maximum (and optimal) = 5. Let's ideal plot of variables is like following: ![enter image description here](http://i.stack.imgur.com/WW2Ya.jpg) I would like to make 5 clusters of this. Something like this: ![enter image description here](http://i.stack.imgur.com/ImEqZ.jpg) Thus I think this is mixture model with 5 clusters. Each clusters have center point and a confidence circle around it. The clusters are not always pretty like this, they look like the following, where sometime two clusters are close together or one or two clusters are completely missing. ![enter image description here](http://i.stack.imgur.com/ingLA.jpg) ![enter image description here](http://i.stack.imgur.com/LJDmZ.jpg) How can fit mixture model and perform classification (clustering) in this situation effectively? Example:               set.seed(1234)     X <- c(rnorm(200, 10, 3), rnorm(200, 25,3),             rnorm(200,35,3), rnorm(200,65, 3), rnorm(200,80,5))     Y <- c(rnorm(1000, 30, 2))     plot(X,Y, ylim = c(10, 60), pch = 19, col = "gray40")


title: Real data collection with small population

text: I performed an actual scientific experiment and I was only able to collect 2 data points. I am wondering what kind of statistical method I can use to analyze the data. I performed a similar experiment before where I was able to collect 4 data points. I calculated the mean, the standard deviation of the population and a confidence interval for this population of N = 4. Could I also calculate the same things for a population of N = 2 ?


title: Whether to stratify or do a simple random sampling from a set of papers to be compared?

text: I have a population/set of papers (~350) that have been categorized into non- mutually exclusive types of in vivo biological target each with different numbers of paper. Another way of categorizing these papers is the method by which some drug against the target is administered. Categories by this categorization are mutually exclusive. I want to compare the merits of the biological target categories. Should I stratify the papers by biological target, or by method of drug administration, or should I just take a simple random sample? Similarly, if I were to do the converse - compare the merits of the methods of administration - how should I sample? Finally, bearing in mind that the method of assessing each category is laborious, what sample size should I be using?


title: Finding the correlation of every row and column between 2 matrices in R, then taking the max and min values

text: I'm using R, and I need to find the correlation between every row and column of matrix A and B (ex: the correlation between the 1st row of matrix A and 1st column of matrix B, 2nd row of matrix A and 1st column of matrix B, 2nd row of matrix A and 2nd column of matrix B, etc.) I realized I can do this by doing cor(matrixA[1,],matrixB) cor(matrixA[1,],matrixB) cor(matrixA[1,],matrixB) ... until I get to the very last row of matrixA (matrixA has 17000 rows) My question is, how do I do this faster without having to type in each command? In addition, the most important part, is how do I get the max and min correlation values of ALL the correlation values I will have calculated in this way? Many thanks! T


title: LibFM sparse data format

text: Does the order of columns in the the sparse format matters in case of libfm ? Can I list the non zero components of X in libfm in any order in a row.


title: What is the reason that a likelihood function is not a pdf?

text: What is the reason that a likelihood function is not a pdf (probability density function)?


title: How to perform a non-equi-spaced histogram in R?

text: From the R docs for hist: > R's default with equi-spaced breaks (also the default) is to plot the counts > in the cells defined by breaks. Thus the height of a rectangle is > proportional to the number of points falling into the cell, as is the area > provided the breaks are equally-spaced. > > The default with non-equi-spaced breaks is to give a plot of area one, in > which the area of the rectangles is the fraction of the data points falling > in the cells. So .. how do I get hist to plot non-equi-spaced breaks? It sounds as if it will calculate the breaks to end up with area one, but I don't see the options. **Edit:** Also, what are recommended ways (in R) to do non-equi-spaced histograms? A typical case would be data that is spiky, causing all the action in one or a few cells, no matter how many are given as "breaks". Another would be two areas of activity separated by a large area of zero, meaning no matter how many breaks, all you see is flat, with two huge narrow spikes. Or perhaps worse, one area of activity, then another area of much less activity far away that causes the graph to be very wide and flat.


title: Model estimation procedure using backward elimination

text: I have run a multiple linear model using `Minitab`. ![enter image description here](http://i.stack.imgur.com/tLGlN.jpg) The result showed that all variables are not significant. So, I use backward elimination. Lastly, I found that one variable is significant, but R square indicates a poor fit. ![enter image description here](http://i.stack.imgur.com/uSRcX.jpg) Can this model be used?


title: Weightining using TraMineR

text: I have read some posts on weighting. However, I am still unclear on the sort of weights I need to use. I am using data from the Longitudinal Survey of Australian Youth (LSAY). This survey provides longitudinal weights for each survey wave (i.e. correction for sampling error and attrition). Because I have weight variables for each period (10) and I can only specify one of them in TraMiner, I am not sure which I should use. I have read that whether weights for the first or last wave should be used. I have not read any reasons why. Can anyone provide me some guidance on this issue?


title: Do only certain pairwise comparisons after significant interaction in two-way ANOVA

text: I am comparing measurements on a test group relative to a control group in three different environmental conditions. I am interested in both differences between environmental conditions and differences between test and control groups. I ran a two-way ANOVA with an interaction term and looked at pairwise comparisons when terms were significant. When the interaction term was significant the Tukey HSD function in R automatically outputs all comparisons. Comparisons between test and test groups on different environmental conditions, comparisons between test and control groups on different environmental conditions and so forth. Needless to say, this resulted in a large number of tests to correct for. My adviser thinks that I should only do three tests to compare test group to control group on each environmental condition (and then only adjust for three tests). I think that because I am interested in differences between environmental conditions in this study, I should run most tests. If I wasn't interested in differences in environmental conditions it should be a nested ANOVA, right? You can see from graphs that the interaction term come from differences between test and control groups on two ecological sites, but it doesn't seem valid to just only run comparisons between groups you 'suspect' will be different. The ones I'm not sure I care about are differences test and control groups on two different environmental conditions. **Is it valid to only run comparisons between groups you are interested in to reduce the number of tests you have to adjust the Tukey HSD p-value for or should you run comparisons on all combinations of groups.** Thanks for help in advance.


title: What sort of learning algorithm to use?

text: I am looking to perform a simple task, > Given restaurant_reviews from week 1 to week n for a particular restaurant,I > look to perform prediction about the number of reviews the restaurant will > get in week n+1.Lets say the number of weeks is X=419 weeks. The problem I > am having is,different restaurants,get reviewed differently, restaurants > that are very popular get 500-600 reviews till week 419. But some > restaurants only get 50-100, hence their reviews are distributed very > differently, I would like to know if an approach of creating competitive > Ensemble Regression Models such that each model is trained on a set of > restaurants of different level of popularity would work for this task? > Example: we might create k bins from our training set based on number of > reviews lets say bin 1 would have restaurants that received between 50-200 > reviews,bin 2 would contain restaurants that received between 201-300 > reviews, bin 3 would contain restaurants that received between 301-400 > reviews and so forth....until bin k. Hence in this case we would get k > regression models each trained on one of the k-bins. I would like to know > whether in this case each of my k models is going to be highly biased or > whether this is an approach even worth trying? Please let me know any better > approaches that I might be able to try. Also any links to papers dealing > with ensemble models for regression would be much appreciated. My apologies for the long post. Thanks in advance.


title: Independent but not identically distributed

text: Let $X_1, X_2,\ldots ,X_n$ be discrete random variables. I'm looking for a way to prove the random variables are independent but not identically distributed. Can anyone suggest some ideas ?


title: Switchpoint detection with probabilistic programming (pymc)

text: I'm currently reading the Probabilistic Programming and Bayesian Methods for Hackers "book". I've read a few chapters and I was thinking on the first Chapter where the first example with pymc consist of detecting a witchpoint in text messages. In that example the random variable to indicate when the switchpoint is happening is indicated with $\tau$. After the MCMC step the posterior distribution of $\tau$ is given: ![enter image description here](http://i.stack.imgur.com/0TGLQ.png) Firstly what can be learned from this graph is that there is a propability of almost 50% that the switchpoint happend on day 45. Though what if there wasn't a switchpoint ? Instead of assuming there is a switchpoint and then trying to find it, I want to detect if there is in fact a switchpoint. The author answers the question "did a switchpoint happen" by "Had no change occurred, or had the change been gradual over time, the posterior distribution of $\tau$ would have been more spread out". But how can you answer this with a propability, for example there is a 90% chance a switchpoint happend, and there is a 50% chance it happend on day 45. Does the model need to be changed ? Or can this be answered with the current model ?


title: Detecting outliers using periodogram?

text: I would like to find outliers in a time series. If there is a an outlier in a time series, how does its periodogram behave? For example I have a time series with 10 elements, and I think 515 in this series is an outlier               475  515 495  507  503   508   496   490   475  510      its periodogram is: ![enter image description here](http://i.stack.imgur.com/cuquc.jpg) if I put 515 on the 7th place:               475   495   507   503   508   496   515   490   475  510      The priodogramm is: ![enter image description here](http://i.stack.imgur.com/Yzs7D.jpg) I have an outlier on different places, but I have two diferent priodogramm, How can I interpret it? or generally   * How can I find an outlier in a time series using priodogramm?   * How can I interpret a periodogram?


title: Random effect with dummy variables

text: By using Stata commands, I tried to run the following three models (5 countries, 3 years):               (1) ols y x1 x2 Dcty1 Dcty2 Dcty3 Dcty4 Dyr1 Dyr2     (2) xtreg y x1 x2, fe     (3) xtreg y x1 x2, re      I would gain the same coefficient of x1 and x2 in (1) and (2) but in (3), it’s naturally different results. Now, I try the following:               (4) xtreg y x1 x2 Dcty1 Dcty2 Dcty3 Dcty4 Dyr1 Dyr2, re      I’m surprised to find the coefficient of x1 and x2 in (4) is the same as in (1). Could somebody tell me why I get these results? In random effect could not put cross-section dummy?


title: Random walk governed by Markov Chain

text: I'm trying to create an underlying markov chain whereby each state governs the probability distribution of a random walk for a stock index. Is there a more formal name for this type of process?


title: what does rating unit means?

text: In a problem in order to calculate the sample size, it is mentioned that: > By use of 7-point Likert-scales, You have identified two segments in the > population that you wish to compare with respect to their mean scores on the > Likert-scale questions. Then it is mentioned some assumptions.i don't understand the meaning of this sentence: > A mean difference of at least 0.75 of a rating unit at the 5% level of > significance. or i can better say, i don't know the meaning of 'rating unit'...


title: ARIMA (0,1,1) or (0,1,0) - or something else?

text: I've just started learning time series so please excuse me if it's painfully obvious; I haven't managed to find the answer elsewhere. I have a data series showing a pretty obvious trend although it's quite noisy. I can take pretty much any division of the data and run classical tests to show a highly significant difference in means. I decided to have a look at time series analysis to see if it could help describe the trend. An ARIMA(0,1,1) model comes out with AIC,BIC=34.3,37.3 (Stata), whilst an ARIMA(0,1,0) model comes out with AIC,BIC=55.1,58.1 - so I understand I'm supposed to prefer the (0,1,1) model. However, the coefficient for the MA(1) is displaying as -0.9999997 (and not showing any p-values). If I try the same in SPSS I get an MA(1) coefficient of 1.000 (I assume SPSS uses opposite signs) with a p-value of 0.990 - does this mean it suggests I drop the term? My understanding is that the effect of a MA(1) coefficient of -1 is basically to remove the old error term and convert the whole series to a linear trend. Does this mean ARIMA is totally unsuitable for my needs? On the plus side it gives me a sensible value for the trend. If I use the (0,1,0) model then I still get a reasonable value for the trend but it's not significant any more. Thanks for your help! EDIT: Thanks for looking in. The trend looks like a fairly linear decrease; the data points seen to fairly noisily rattle around above and below a trend line. The ARIMA (0,1,1) model produces something that's not far off a straight line decrease which seems sensible - the (0,1,1) produces what is essentially a lagged version of the data, translated down by one month of trend. The data aren't stationary (due to the trend) - though the first differences seem to be. I don't think the (0,1,1) is a bad model - I'm just a little confused by the p-value seeming to suggest I should drop the MA term - or wondering if it means I should bin ARIMA entirely! EDIT2 @vinux - thanks for the suggestion; that makes a lot of sense (and seems to be what the -1 MA term is trying to create?). I've uploaded as many graphs as I could think of as people had requested. ![tsline y - graph of raw values](http://i.stack.imgur.com/OR2sj.png) ![tsline D.y - graph of differences](http://i.stack.imgur.com/xlzl2.png) ![ac y - autocorrelations of y](http://i.stack.imgur.com/8lQ57.png) ![pac y - partial autocorrelations of y](http://i.stack.imgur.com/IlqiP.png) ![ac D.y - autocorrelations of first differences](http://i.stack.imgur.com/gUa6p.png) ![pac D.y - partial autocorrelations of first differences](http://i.stack.imgur.com/KE5Sa.png) I've also put the monthly data up in CSV format at pastebin


title: Why do I not get a p-value from this ANOVA in R?

text: Here's the data:               > tires <- data.frame(Wear  = c(17, 14, 12, 13, 14, 14, 12, 11,                                     13, 13, 10, 11, 13, 8, 9, 9),                           Brand = rep(LETTERS[1:4], 4),                           Car   = as.character(as.roman(rep(1:4, each = 4))))     > tires        Wear Brand Car     1    17     A   I     2    14     B   I     3    12     C   I     4    13     D   I     5    14     A  II     6    14     B  II     7    12     C  II     8    11     D  II     9    13     A III     10   13     B III     11   10     C III     12   11     D III     13   13     A  IV     14    8     B  IV     15    9     C  IV     16    9     D  IV      Now I fit a two-way ANOVA with interaction:               two.way <- aov(Wear ~ Brand + Car + Brand:Car, data = tires)      Finally, no p-values:               > summary(two.way)                 Df Sum Sq Mean Sq     Brand        3  30.69  10.229     Car          3  38.69  12.896     Brand:Car    9  11.56   1.285      A regular two-way ANOVA (i.e., `Wear ~ Brand + Car`) gives me p-values:               > summary(aov(Wear ~ Brand + Car, data = tires))                 Df Sum Sq Mean Sq F value  Pr(>F)        Brand        3  30.69  10.229   7.962 0.00668 **     Car          3  38.69  12.896  10.038 0.00313 **     Residuals    9  11.56   1.285                        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Is there a way to interpret this? The interaction plot shows me that there is definitely interaction between `Brand` and `Car` so I am hoping to incorporate this into my model.


title: Show Regression with Arima Errors Equivalent Form of Differenced Variables

text: How can you show that the regression $y_{t}=\beta_0 + \beta_1x_{t}+\eta_t$ where $\eta_t$ is arima(1,1,1) is equivalent to $y'_t = \beta_1x'_{t}+\eta'_t$ where $\eta'_t=\phi_1\eta'_{t-1}+e_t$ and $'$ denotes a first difference? Also, is it possible to construct a model where everything is differenced except $\eta_t$, as in $y'_t = \beta_1x'_{t}+\eta_t$ where $\eta_t=\phi_1\eta_{t-1}+e_t$ where $'$ denotes a first difference or does the error term always have to be differenced if y or x is?


title: Mixing two linear regression models

text: In time series analysis, I have one predictor $X_1$ that has a higher $R^2$ when regressed against $Y$ sampled at **10 minutes** interval. Another predictor $X_2$ has a lower error when fitted against $Y$ at **20 minutes** intervals. The predictors are built from different data sets and unrelated to each other. For example, assuming we run two different regression models for $X_1$ (10 minutes horizon) and $X_2$ (20 minutes): at time $t_0$ the $X_1$ model predicts 10 and $X_2$ 30. Ten minutes later at time $t_1$ model $X_1$ forecasts 100. So now, at time $t_1$, we have $X_1$ predicting 100 and the **previous** $X_2$ forecast calculated at $t_0$ predicting 30 effectively for the same horizon target $t_2$ 10 minutes out. How do I address the two models overlapping horizons? I don't necessarily need to build two models, I'm trying to understand how this kind of time series problems can be dealt with.


title: Proc Mixed statement

text: Suppose I have a data set:               Geography Quarter  Y X1 X2     A           1      3  1  4     A           2      2  4  4     A           3      5  5  2     A           4      1  3  4     B           1      3  1  6     B           2      8  6  7     B           3      2  4  1     B           4      9  7  3      What does this statement do?:               proc mixed data = data;     model Y = X1 X2;     run;      Does it fit one model or multiple models? I am trying to understand what proc mixed actually does with this dataset. Would this statement simply fit one model for all data points? And would adding a class statement create a model for A and a model for B?


title: Conjugate of Weibull with shape known

text: This isn't exactly a homework problem but rather a self-selected problem I'm doing to prepare for a midterm. I can see from Wikipedia that it is an inverse gamma but I am unable to reach the posterior suggested by the parameter updates. The Weibull as given on Wikipedia is: $f(x_1,..,x_n|\lambda, k)=\frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}\text{exp}\left\\{-\left(\frac{x}{\lambda}\right)^k\right\\}$ Then $L(\lambda|x_1,\dots, x_n)\propto \prod_x \lambda^{-k}\text{exp}\left\\{-\frac{x^k}{\lambda^{k}}\right\\}\\\ =\lambda^{-nk}\text{exp}\left\\{-\frac{\sum_x {x^k}}{\lambda^{k}}\right\\}$ The inverse gamma prior is $P(\lambda|\alpha,\beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}\lambda^{-\alpha-1}\text{exp}\left\\{-\frac{\beta}{\lambda}\right\\}$ The posterior then is $P(\lambda|x_1,\dots,x_n)\propto L(\lambda|x_1,\dots, x_n) P(\lambda|\alpha,\beta)\\\ =\lambda^{-nk}\text{exp}\left\\{-\sum_x {x^k}\lambda^{-k}\right\\}\lambda^{-\alpha-1}\text{exp}\left\\{-\frac{\beta}{\lambda}\right\\}\\\ =\lambda^{-nk-\alpha-1}\text{exp}\left\\{-\frac{\sum_x {x^k}}{\lambda^{k}}-\frac{\beta}{\lambda}\right\\}$ The problem that I am left with is to coerce an update for the $\beta$ hyperparameter that does not contain $\lambda$. The $\alpha$ also looks different from Wikipedia but I am going with explaining that as the difference between the posterior with observations and without, though honestly I'm confused about that as well. Wikipedia simply has the summation, $\sum_x {x^k}$ as the update for $\beta$ and $n$ for the update of $\alpha$ (not $nk$). All of my efforts so far to coerce that $\beta$ update have failed.


title: Outlet for short communication

text: I have written an 8-pages technical paper (more theoretical than applied) in statistics. I only know of two journals (Biometrika and Statistics and Probability Letters) that publish such short communications and I was wondering whether there are other journals that take this kind of articles.


title: Adding the variances of spatially correlated data

text: I have a geographical grid, for every cell in that grid I know that the variable I am looking at is normally distributed and I have its mean and variance. Now, imagine I am looking at an individual grid cell for which I have little information, thus I would like to infer the standard deviation for that grid cell from the 8 surrounding cells. I know variances can be summed for independent variables, which is obviously not the case here because they will be spatially correlated. I have little knowledge of geostatistics, although I am aware they could be relevant for this. I am happy to look in that direction as well but not sure where to start. Any pointers would be great.


title: Inconsistent x-y relationship from paired t-test and linear regression

text: I am analyzing an experiment comparing the effect of treatment A vs. B on the matched subject. Here are the measurements on 34 subjects:                A     B     -1.15 -1.16     -1.13 -0.94     -0.16 -1.18     -0.37 -1.20     -1.09 -1.20     -1.20 -1.20     -0.94 -1.20     -0.84 -1.16     -1.18 -1.17     -1.20 -1.11     -0.78 -0.68     -0.83 -0.73     -1.05 -1.20     -0.71 -1.20      0.07  0.12     -1.20 -0.98     -1.20 -1.20     -1.02 -1.17     -0.28 -0.84      1.33  1.47     -1.19 -1.20     -1.20 -1.17     -0.40 -1.20      0.66 -0.21     -0.63  0.21     -0.88 -1.16     -0.46 -1.20     -0.76 -1.20     -0.38 -1.20     -0.67 -0.97     -0.90 -1.20     -0.90 -1.20     -1.20 -1.15     -1.01 -0.79      The differences between the two treatment (dat[,"A"]-dat[,"B"]) looks normally distributed. I first applied a paired t-test:               t.test(dat$A, dat$B, alternative = c("two.sided"), mu = 0, paired = TRUE)          Paired t-test     data:  dat$A and dat$B      t = 2.894, df = 33, p-value = 0.006692     alternative hypothesis: true difference in means is not equal to 0      95 percent confidence interval:      0.05870022 0.33659390      sample estimates:     mean of the differences                    0.1976471       The paired t-test indicates that on average treatment A has a significantly higher measurement than treatment B. On the other hand, I applied a linear model on A~B:               mod1 <- lm(A ~ B, data=dat)          Coefficients:             Estimate Std. Error t value Pr(>|t|)         (Intercept) -0.03345    0.12739  -0.263    0.795         B            0.75111    0.11801   6.365 3.79e-07 ***      The 95% confidence interval of the coefficient for B (0.51-0.99) does not cover 1. This result indicates that treatment A has on average a smaller measurement than treatment B, which is contradictory to the findings from the t-test. Can anyone help me to explain these contradictory findings? and to expand my question: Does a paired t-test (test mean of difference against 0) equal to a linear regression without intercept (test the coefficient against 1)? I mean in terms of testing against the null hypothesis, rather than the estimate or the meaning of the coefficient. Because both tests are testing against the null hypothesis that $B_{i}-A_{i}=\epsilon_{i}$.


title: How to compare the outputs of two algorithms computing SVD?

text: For example we have 2 algorithms from R: SVD and irlba and I want to compare them int terms of speed,memory and precision. But I don't understand how to compare output of algorithms, they must be differ on some epsilon, but how to define epsilon? And why sign of principal vectors can be different? Here is some test code:               rows= 10     cols= 8     A <- matrix(runif(rows*cols),rows,cols)          pc=5     a <- irlba(A,nu=pc,nv=pc)     b <- svd(A,nu=pc,nv=pc)          #false     a$u == b$u     a$v == b$v          #show all vars     a$u         b$u     a$v         b$v     a$d         b$d          #compare u     .list <- list(abs(a$u), abs(b$u))     Reduce('-', .list)          #compare v     .list <- list(abs(a$v), abs(b$v))     Reduce('-', .list)          #compare d     .list <- list(abs(a$d), abs(b$d[c(1:pc)]))     Reduce('-', .list)      output:               > #compare u     > .list <- list(abs(a$u), abs(b$u))     > Reduce('-', .list)                    [,1]          [,2]          [,3]          [,4]         [,5]      [1,]  6.383782e-16 -5.551115e-17 -1.082467e-15  1.665335e-15 -0.123009729      [2,]  0.000000e+00  5.689893e-16 -3.053113e-16 -1.665335e-15 -0.117362383      [3,]  1.110223e-16 -5.481726e-16 -2.220446e-16  7.216450e-16  0.142286325      [4,]  3.885781e-16 -3.053113e-16 -1.526557e-16 -2.942091e-15  0.096784035      [5,]  4.996004e-16 -9.436896e-16  1.776357e-15 -1.956768e-15 -0.081504771      [6,] -2.775558e-16 -6.661338e-16  7.216450e-16  1.137979e-15  0.007450404      [7,] -5.551115e-16  4.440892e-16 -5.273559e-16  5.551115e-17  0.018699113      [8,] -5.551115e-16  7.216450e-16  4.440892e-16  1.665335e-16 -0.013035052      [9,] -4.440892e-16  1.054712e-15 -1.942890e-16  2.498002e-16  0.048095160     [10,]  3.330669e-16  3.330669e-16  7.771561e-16  6.106227e-16  0.044937739     >      > #compare v     > .list <- list(abs(a$v), abs(b$v))     > Reduce('-', .list)                   [,1]          [,2]          [,3]          [,4]        [,5]     [1,] -1.887379e-15 -2.775558e-16 -9.020562e-17 -1.193490e-15 -0.02102106     [2,]  1.831868e-15 -1.276756e-15 -6.661338e-16  3.386180e-15 -0.09872305     [3,]  7.216450e-16 -3.885781e-15 -9.436896e-16  9.992007e-16  0.04699264     [4,] -2.053913e-15  2.331468e-15 -4.996004e-16  1.498801e-15  0.05876057     [5,]  8.881784e-16 -5.551115e-16 -8.881784e-16 -1.942890e-16  0.16146650     [6,] -2.053913e-15  8.881784e-16  5.551115e-16 -2.775558e-16  0.11420391     [7,]  1.498801e-15 -1.665335e-16  1.665335e-16  9.714451e-16 -0.13600871     [8,]  8.326673e-16 -2.831069e-15 -7.042977e-16 -2.997602e-15 -0.01507590     >      > #compare d     > .list <- list(abs(a$d), abs(b$d[c(1:pc)]))     > Reduce('-', .list)     [1] -8.881784e-16  0.000000e+00  8.881784e-16  2.220446e-16 -2.483630e-02      here is also some comparision https://github.com/graphlab- code/graphlab/issues/48


title: Monte Carlo Tree Search : Expansion Step

text: The algorithm asks to **"Start at root node, recursively select optimal child nodes until a leaf node L is reached."** , and then **"If L is a not a terminal node (i.e. it does not end the game) then create one or more child nodes and select one C."** My question is which node should we expand? And again, suppose the leaf L has children _C1, C2, ..., Cn_ and we expand _C1_ in this step. So L is not a leaf anymore. Then in future how we can expand other _Ci_ 's when the algorithm asks to expand when we hit a leaf? Please someone explain the **expansion** step.


title: Naive vs. robust variance in GEE

text: I'm analyzing a data set including 40 participants. Each participant played a decision-making game 80 times. I am using a GEE analysis clustered by participant, with an exchangeable working matrix. What I am unsure about is whether I should use the naive S.E (in Stata this is the default option), or work with robust S.E. I am under the impression that the robust option is "too harsh" on the data. I would happy to hear your thoughts on the subject, and as always, any reference is welcome. Thanks!


title: Statistics Test for bivariate frequency table

text: First time user here! I have some troubles in finding some statistic method that would fit my purpose for testing bivariate frequency table. Below is a sample frequency table I have created with ordinal data. ![enter image description here](http://i.stack.imgur.com/CffKV.png) The rows and columns are labeled with ranks from 1 to 4. As you can see, most of my data is concentrated at the bottom-left corner for the sample. My goal is to quantitatively define this characteristic and perhaps perform a statistical test for estimating the dispersion in the population. I have tried to the Chi-squared contingency table and Fisher's exact, but my table contains a few zeros cells. Please advice any other appropriate test that I should attempt. Thank you. **ADDed: The experiment is designed as follow,** Four stocks are picked and their expected return and standard deviation are measured for 40 observations/periods. Within each observation, rankings are assigned to each stock for both y (expected return) and x (standard deviation) as 1 being the lowest and 4 being the highest. The frequency table above focuses on one of the four stocks. if most of the observations fall into region of (y >= 3 & x <=2), I would characterize this stock as "low risk, high reward" option.


title: What are finite window effects?

text: I'm reading a paper that uses a Poisson process to model real world events. The authors mention "finite window effects". What are finite window effects? Here is quote from the paper where the authors first mention the term: > If the data come from a Poisson process, then a histogram of inter-event > times will be roughly uniform when looking at a short inter-event time > window due to minimal ﬁnite window eﬀects. This is because, with an inﬁnite > window, inter-event times from a Poisson process are uniformly distributed.


title: Absent categorical data levels in Bootstrap samples

text: I have a huge dataset ($n$ around five million, $p$ around three thousand) for a classification problem, where my interest is predictive class probabilities for test data, not the target. I shall be using bootstrap samples to estimate the smoothed probabilities from various models. The problem is when I bootstrap the cases not all the levels in the categorical predictors will be in every resample. Given that I need parameters for all levels when working with the test set, I'm unsure of what to do. I read Agresti's _Categorical Data Analysis_ , but there doesn't seem to be mention of this. I have thought of 2 possibilities:   1. Insert a base composition of varied cases in every resample so that all levels are included for all predictors.   2. Define all levels for each categorical variable with reference to the data schema and run the models. I'm having to think about this because my resamples are pulled as CSV data from a SQL database and I would usually use `read.csv()` which automagically handles levels, preference relations and levels for categorical data using the cases in the CSV file. (Pulling all of the data in one csv is not an option due to resource constraints.)


title: PDFs and probability in naive Bayes classification

text: I have seen a few times the technique of using the Gaussian PDF for continuous features in Naive Bayes. here and here. Illustrated in the first link: ![enter image description here](http://i.stack.imgur.com/NbRty.png) How is this possible? I always learnt that the PDF is not a probability -- as the probability of any exact value of x is zero.


title: SVM and non-linear predictive models - feature selection

text: Just throwing out a general question. What do people think of applying feature selection methods when using SVMs to build predictive models? I understand that SVM have built in regularization with how they're trained, but I've heard for certain cases (e.g. features >> examples) feature selection is still useful. I'm mainly talking about non-linear SVMs (e.g. kernelized). Is there a difference in effectiveness of feature selection when say # features >> # examples vs. # features ~ # examples vs. # features << examples? If you don't think feature selection is useful, please explain why. If you think feature seleciton is useful (and have practical experience in this regard), please let us know what methods you've used and how they worked. Also open to entirely different approaches to non-linear predictive modeling other than SVMs. Thanks!


title: iOS App Store Data

text: I'm looking for some data sources for apps and/or transactions in the iOS App Store so I can perform some analytics. Any database with a dump of titles, prices, description text, reviews, etc. would be amazing. There are some out there that I have found but it is difficult to figure out how credible they are. Thanks!


title: Odds ratio and confidence interval in meta-analysis

text: I need to calculate the pooled odds ratio and associated 95% confidence interval for meta-analysis of 2 studies about the risk of bleeding. The only information I have is the odds ratios and 95% confidence intervals. They are 2.7 (1.8 – 4.0) in the first, and 1.3 (0.5 – 3.4) in the second study. I computed the standard errors, weights, and pooled SE, OR and CI from the available ORs and CIs. According to my calculations, the standard errors are 0.204 in the first, and 0.489 in the second study, which gives the pooled OR and CI of 2.49 (1.72 – 3.60). I’m not sure about this and I would appreciate if someone checks this up. Thanks in advance!


title: How many times can I guess a number from Sample 1 belongs to Sample 2 before concluding that they are not the same?

text: First of all, I apologize for my half-forgotten understanding of statistics. I hope I can accurately describe what kind of test I'm after in a way that doesn't sound confused. I have one sample and one population. I would like to know the probability that my sample is part of that population. But I also want to know what are the maximum number of "guesses" from my sample it will take in order to conclude that it is, or is not a part of that population. Let's say I'm looking for a 5% confidence that the sample is or is not part of the population. So basically, I take a number randomly from my sample, and check to see if it exists in the population. How many times should I guess until I can feel confident that the sample is or is not part of the population? I have a funny feeling this is a pretty elementary question and I'm just confused by the word problem aspect of it. Anyway, I know the size of the sample and the size of the population, if that helps. A mean of the sample & population is meaningless because these numbers are values, they're labels. I'm certain there's a more correct way to describe this. But essentially the numbers are like house addresses, so the difference between 1 and 5 compared to 25 and 300 is meaningless. Thanks for any help. I suspect this problem is related to statistical process control. Thanks again.


title: Are the law of statistical regularity and inertia of large number applicable at level of meta- analysis?

text: Principles of sampling are applied on study of sample data on variables collected from number of observations(n) for estimation of population parameter(effect size). Meta analysis is applied on sample correlations(r's)/effect sizes(d's) of number of studies(k).


title: The mean and Variance of $log(Combination Index)$

text: In biological sciences, ($CI$) stands for the Combination Index which is a conventional method for dose-response assessment and drug interaction analysis. It can be defined as: $$CI = \frac{\bar{D}_{1} \hat{\beta}_{1}}{ln(\bar{SF}_{1}) - \hat\alpha_{1}} + \frac{\bar{SF}_{1} - \hat{\alpha}_{2}}{{0.5} - \hat\alpha_{2}}$$ Notes:   1. $\bar{D}_{1}$ is independent of $(\hat\alpha_{1}, \hat\beta_{1})$.   2. $(\hat\alpha_{1}, \hat\beta_{1})$ are independent of $(\hat\alpha_{2})$.   3. $log(SF_{1}) = \alpha_{1} + \beta_{1} x + \epsilon$   4. $E(\bar{D}_{1}) = \mu_{D}$, and $Var(\bar{D}_{1}) = \sigma^{2}_{D}$ I'm trying to find $E({log(CI))}$ and $Var(log(CI))$. I know this can be done using delta method, but I'm not sure how... your help is appreciated. Thanks!


title: Define a hypersherical neighborhood around a target point with radius r'

text: _This is a homework problem that I need some help with._ > QUESTION: We assume that N data points are uniformly distributed in a > 100-dimensional unit hypersphere (i.e. $r = 1$) centered at the origin, and > the target point $x$ is also located at the origin. Dene a hyperspherical > neighborhood around the target point with radius $r_0$. This is how I understood the problem: ![drawing of the problem](http://i.stack.imgur.com/ol5V6.jpg) So my question is: **to find the neighborhood with radius $r'$, do I have to compute the median distance of the closest point from the origin, which is then given by $$ {\rm distance}(d,N) = \bigg(1-\left( \frac{1}{2} \right)^{\frac{1}{N}}\bigg)^{\frac{1}{d}}\ , $$ or is it given by the surface area of the hypersphere with radius $r'$?


title: Obtaining a Probability Distribution From a Survival Function

text: **Edit** : I basically want to have a probability curve where a X value of 0.002 would be associated with a Probability of 1 and would also have data points of (0.005,0.1), (0.008,0) which is seen in the graph below. However, I then realised I would need to check the probability that the aforementioned random variable is greater than another random variable, which is normally distributed, and for that I would need both to be distributions. So, is there some way I could create a PDF having a CDF like the graph below? I realise the CDF must always be increasing, so would it be possible to have a graph increasing with the x values increasing from 0.008 to 0.002? **Edit** : ![enter image description here](http://i.stack.imgur.com/z57Tj.jpg)


title: Sample data from combination of two probability distributions

text: I want to generate a mock catalogue. I have access to two independent sets of real data and I want to use their properties to make the mock catalogue: The first one contains the information from **magnitude** and **redshift** ($z$) of galaxies, while the second data set has **magnitude** and **weighted ellipticity** of galaxies. In my final simulated catalogue, I want to generate **redshift** and **weighted ellipticity** for a set of randomly generated galaxies. I would like to use the underlying properties of **magnitude** versus **redshift** and **magnitude** versus **weighted ellipticity**. However I also have a well- define formula that discribes the distribution of redshift and the functionality is given as following (`GeneralRandom` is a class that is able to generate a random sample follows the given distribution $p(z)$):               def p(z):         z0=1./3.;eta=1.0         value=eta*(z**2)*np.exp(-1*(z/z0)**eta)/scipy.special.gamma(3./eta)/z0**3          return value          catalogue_generate=GeneralRandom(x=np.arange(0.0, 1.5, .001), p_func=p, Nrl=10000)      ng=24000     catalogue_generate.set_pdf(np.arange(0.0, 1.5, .001), catalogue_generate.p_val, 10000)     redshift=catalogue_generate.random( ng)[0]     redshift=random.sample(redshift, ng)      The sample with input **redshift distribution** is produced by above procedure. I also used `scipy.stats.gaussian_kde` to estimate the probability distribution of **magnitude** versus **redshift** and **magnitude** versus **weighted ellipticity** for two available data sets. Now I am confused how I could combine the `kde` computed from two distribution and obtain just my interested parameter in the mock catalogue, `weighted ellipticity` from the `sampled redshift`? The plot is illustrating the Gaussian KDE estimated for two data sets. I don't know how I could independently use the output of `stats.gaussian_kde` to just insert computed `redshift` and obtain `magnitude` and then from estimated `magnitude` before, and estimated second KDE estimate the `weighted ellipticities`?![enter image description here](http://i.stack.imgur.com/aIoBQ.jpg) Is it statistically correct to say the problem looks like this: $$P(w_\epsilon)=\int P(Mag|w_\epsilon)\int P(z|Mag)P(z)dz$$ I also don't know computationally, how it should be done?!! I will appreciate if somebody can describe how I could tackle my problem and a bit of python programming tips or introducing useful libraries.


title: Does this data type require this manual standardisation, or can this be circumvented using a random forest classifier?

text: I have some data that represents students exam scores, this data is used to predict what rank the student will be in 'todays' exam. Some of the data I use in order to predict this is the previous exam scores that these students achieved. Unfortunately the previous exam scores are not always from the same subject as todays exam. For instance todays test may be calculus and yesterdays (or the last one taken) may have been Latin. If I can describe my previous standardisation methodology to help add some context. First thing to note is that I am not currently trying to work out what test score a given pupil should achieve for this exam they are taking, I am trying to calculate where they should rank against their peers within that exam class.   * My original method: I calculated what the average score in a calculus exam was by a student that achieved a 98 in the Latin exam (which he took as his last exam) via a simple non linear regression. I do this for all the students sitting a particular exam, calculating what I think their previous exam score would be worth if it was in todays subject. I even do this for same subject scores, the reason being when a student scores a ridiculously high score one exam, in reality on average they achieve slightly lower in the next exam (according to my regression.)   * What I then do is find the difference between the largest last time exam score for a class and the last time exam score of every pupil taking todays exam. e.g. The highest scoring pupil will have a value of `0` and those who achieved worse will have a positive number reflecting how far away they were from the person who achieved the best score, lets call this value `deltaE` for reference purposes.   * Then I calculating subsequent percentage exam class win of different `deltaE` values, for example those with a `deltaE` of 0 are found to have a 18% chance of ranking highest in their next exam.   * I then use this data in an attempt to predict where a student will rank using a random forest classifier, including all three parameters 18,100,0.18 representing amount of times a `deltaE` of 0 rank top in their next exam, from a total number of `deltaE-0` students and finally the 18/100 represented as a decimal This method of standardising the different subjects into an expected subject relevant score to me seems messy and more than likely the long way of doing things. I am wondering if it is possible for me to use the random forests functionality to remove the need for this standardisation? For instance could I use a numerical representation of each subject and include this in my data? For instance Latin could be 18,100,0.18,2 however if the subject was calculus it would be 18,100,0.18,1 or something to this effect. All comments and suggestions welcome, someone with experience of random forests would be ideal to comment on this, but if your experience is in other modelling programming backgrounds feel free to put your thoughts across. Thanks in advance


title: Basics of data analysis with exponential data/noise?

text: Can you suggest a tutorial or book chapter on basics of data analysis with exponentially distributed data / exponential noise, at undergraduate level?   By "basics" I mean:   * is a given set of data exponentially distributed?   * should one generally use median instead of mean for exponential data? Should one simply trim the top 10% (for some value of 10)?   * least-squares curve fitting, where I believe errors to be exponential? (Ok, I'll be made to write out "it's not that simple" 100 times.)


title: How to analyse three categorical variables

text: I need some help identifying a test to use for three categorical variables: Subject (maths, business etc), Big 5, and Learning style. I am carrying out research on whether there is a relationship among the above three variables. There are no scores, only categories. I looked at the chi-squared test, but it doesn't appear to be helpful if I have more than two variables. Participants are all one gender.


title: Standard error help answer

text: The standard population error of the sampling distribution when we know the population standard deviation is equal to $σ_M = σ/\sqrt n$, or is it, $S_m = σ/\sqrt n$?


title: Confidence interval for the standard deviation on a bimodal distribution

text: I have a bimodal distribution, and I wish to estimate the mean and standard deviation of the population (well, these could be 2 sub-populations according to the shape). With the mean I have no problem, the sample is relatively large (n~70), and so the sample mean is approximately normally distributed according to the central limit theorem. I can calculate the corresponding 95% CI. Now I want to do the same for the standard deviation, but I am not sure if the same theory apply. Can I say that the variance is approximately chi squared distributed when the sample is large ? Edit: To test my question, I ran this simulation in R:               mu1 <- 23        mu2 <- 40     sig1 <- 4.3     sig2 <- 2.8     cpct <- 0.32     m <- 10000     means <-rep(0,m)     sds <-rep(0,m)     vars <-rep(0,m)          bimodalDistFunc <- function (n,cpct, mu1, mu2, sig1, sig2) {       y0 <- rnorm(n,mean=mu1, sd = sig1)       y1 <- rnorm(n,mean=mu2, sd = sig2)            flag <- rbinom(n,size=1,prob=cpct)       y <- y0*(1 - flag) + y1*flag      }               for (i in 1:m)     {         bimodalData <- bimodalDistFunc(n=10000,cpct,mu1,mu2, sig1,sig2)         means[i] = mean(bimodalData)             sds[i] = sqrt(var(bimodalData))         vars[i] = var(bimodalData)     }          par(mfrow=c(2,2))     hist(bimodalData)     hist(means)     hist(sds)     hist(vars)      And I got the following result: ![enter image description here](http://i.stack.imgur.com/AZq6U.jpg) The distribution of the variances is not chi squared, but why is it normal ? (the bimodal histogram is an example of one of the simulated distributions)


title: How would you approach weighting of open response lists in a statistical manner?

text: How would you deal with lists of varying lengths? For example, if you asked for an _ordered list_ of your favourite movies without any constraints how would you value the weight of someone who submits 5 samples and another who submits 20 samples? Aside from the linear valuations, crudely you could give each participant's first selection a score of n, 2nd = n-1, 3rd = n-2... nth=1/nth value (n being the longest list supplied) and build an aggregated score across the sampled population. But that feels far too crude for proper analysis. Is there a better weighting distribution with a tapered decay to suit a long tail? Instinctively, the person who supplies a longer submission should weigh higher than a shallow submission. But I'm unsure of the social research techniques that might be suitable. Is there a statistical model that addresses this kind of analysis for open response lists?


title: testing independent of data assumption, before doing Poisson Regression

text: I have got a set of data and decided to do a Poisson regression so $log(\lambda)=X\beta$. I wonder if there is any test available for such purpose? The first thought is to use the Durbin-Watson test statistics to test the independence of data. My primary concern is that since the mean and variance is not the same for each observation so if the Durbin-Watson statistic is correct? If not, is there any suggestions? Thank you!


title: Testing presence of one outcome over another

text: I have what I think is a simple problem, but I can't seem to think of how to test it. I have a sample of 2649 patients, of whom 125 have a specific disease. The disease is classified into "early" and "late", and I want to test the difference between the probability of being "early" vs "late", and, if possible, calculate a relative risk. I can get the probabilities easily enough, and I can get the ratio of the probabilities to get a relative risk. My question is how do I test these variables or get confidence intervals: since both outcomes are from the same population the equations I'm used to don't work. Thanks in advance, I've included some simple sample R code below.               sample = c(rep(0,2517),rep('A',68),rep('B',64))     p = c(sum(sample=='A'),sum(sample=='B'))/length(sample)     RR = p[1]/p[2]


title: Adding errors in quadrature for mean values

text: I understand that for $ Z = \frac{X}{Y}$, the error in quadrature would be calculated as, $\Delta Z = Z \sqrt{\big( \frac{\Delta X}{X} \big)^2 + \big( \frac{\Delta Y}{Y} \big)^2}$. In my case, X and Y are the mean intensities of some regions of interest (ROI) within an image. Therefore, would the $\big( \frac{\Delta X}{X} \big)$ correspond to the standard deviation of the voxels divided by the mean of the voxels within the ROI?


title: How to convert a dendrogram back into a distance matrix?

text: Example code:               our_dist <- dist(USArrests[1:4,])     dend <- as.dendrogram(hclust(our_dist , "ave"))     plot(dend )      I would now wish to have a "dend2dist" function which will turn dend back to our_dist. Of course, it is not possible to do in full (AFAIK), since the process of turning the dist into the dendrogram aggregated a lot of the data which can not be retrieved. But still, I would like for our "best guess" (assuming which aggregation method was used for the dendrogram), as to what is some distance which would be able to reproduce the original dendrogram. And of course, tips on how to best implement such a function in R would also be nice. Any suggestions will be most welcomed, thanks.


title: Variety of ways of reporting aged debt

text: I have been looking a different way to do reporting on Aging report customer Account balances (A/R balance). The A/R balance is presented in terms of age of debt to measure the total open debt-by-debt ages in the A/R system. The age of a debt is determined by its due date, and the number of days past the due date, for example 1-30 days, 31-60 days, and so forth. Currently this is what I been graphing as an example described above 1-30 days, 31-60 and last bucket is 120+ and have found it boring and wanting if there are other ways doing this reporting. Here are a few data points:               day.past.due    amt      57              46,207.74      321             25,275.76      625             21,043.17      2618            26.00      650             2,952.40      465             484.22      1407            10,090.55      979             239.10      235             25,965.41      218             50,715.57      366             15,616.98      608             15,819.15      534             4,087.55      695             7,138.07      844             3,680.72      684             24,511.03      1797            1,786.16      985             7,004.80      108             38,570.32      258             22,897.39      148             24,200.36      461             13,240.68      444             9,386.73      489             2,665.33


title: C5.0, boosting, and mislabeled data

text: I'm trying to model a binary classification problem. 5 continuous features, slightly imbalanced dataset (about 60 in one class, a little more than 200 in another). So far I've tried kNN, LDA, and C5.0. C5.0 works very well (e.g. the AUC for ROC is very good). However, I've read that boosting a la AdaBoost has problems in the presence of mislabeled data, and that C5.0's boosting is reminiscent of AdaBoost. My data is mislabeled insofar as perhaps 5% of them are probably wrong. (I'm vague here, because it's a real-world problem where the ground truth is fuzzy--- good, but not perfect.) Does anyone know how well C5.0 boosting performs in the presence of mislabeled data?


title: Main idea of Bagging

text: I just read this post and several other websites, but I still don't understand what bagging is. I understand it is an algorithm for machine learning, that it improves stability and accuracy of the algorithm and decreased the variance of my prediction, but what is the main idea behind this algorithm? Is it to put a dataset into a bag? _For example, to me the main idea behind boosting is to boost records that are weighted incorrectly._


title: How to decompose a model using linear regression?

text: I'm using "decompose" in a figurative sense here. I have a simple regression model: y = b + b1x1 + b2x2 + e The "e" term is the residual and "b" the intercept. I want to show the contribution of x1 (and also x2) toward y For the contribution of x1 I remove x2 by setting it to zero. That yields: y|x1 = b + b1x1 For x2 then: y|x2 = b + b2x2 The problem is that these two contributions do not sum to the modeled y values as the intercept is double-counted: y|x1 + y|x2 = 2b + b1x1 + b2x2 != y When I force the regression thru the origin it works fine and there are solid physical reasons to do so but I am looking for a more elegant way to do this? Comments to answers here (I can't add comments below for some reason): To Peter, yes, that is it. Imagine a stacked bar chart. I want my bar chart to have 3 components, the amount of y determined by x1, the amount of y determined by x2, and the residual. These should add up to the raw data.


title: What are the differences among latent semantic analysis (LSA), latent semantic indexing (LSI), and singular value decomposition (SVD)?

text: These terms get thrown around together a lot, but I would like to know what you think the differences are, if any. Thanks


title: Estimate single ARIMA for multiple timeseries

text: I have two groups of time-series, each group represents one type of data. However within each group, each time series may be fitted with a different ARIMA(p,d,q) from the other time series in the same group. I need to create a single model for each group (`Model_group1`, `Model_group2`). I tried the approach mentioned by Rob Hyndman in: Estimating same model over multiple time series. I need to use these two models to classify any time series to one of these two groups. For each time series, I calculated the AIC of `Model_group1` and `Model_group2`, and the model with smaller AIC will mean that the time series belongs to its corresponding group. I have three problems:   1. I received a warning message               Series: ts      ARIMA(3,0,2) with non-zero mean           Coefficients:              ar1     ar2     ar3     ma1      ma2  intercept           0.0714  0.1417  0.0000  0.0893  -0.0871     0.1169     s.e.     NaN  0.1381  0.0127     NaN   0.1436     0.0026          sigma^2 estimated as 0.2202:  log likelihood=-33822.63     AIC=67659.26   AICc=67659.26   BIC=67725.99     Warning message:     In sqrt(diag(x$var.coef)) : NaNs produced      This message was returned by only one of the group models. Does that mean that the fitted model is not correct?   2. I got two different results using               auto.arima(ts, allowdrift=FALSE, stepwise=FALSE)     auto.arima(ts, allowdrift=FALSE, stepwise=TRUE)        3. When I tested the resulting models, the majority of the time-series were classified as `group_1`, even when I test one of the time series used to build the long time series of `group_2`. I need to mention here that the composed time series of `group_1` is quite shorter than the time series of `group_2`. Are there any expected reasons for that?


title: Comparing patient proportions between one ER and the national average

text: I am looking to compare the proportion of emergency room attendances by hour of the day against a national average proportion (plus some other tests for ER variables such as referral methods and discharge destinations). Ho = no difference between the two distributions.   I was intending to use a paired t-test for the normally distributed data and a Wilcoxon signed rank test for the non-normal. I have never used the Wilcoxon using proportions, though. Is this acceptable and if so does anyone know a good text / previous research to read? Many thanks.


title: Canonical Correlation analysis without raw data

text: I want to run a Canonical Correlation (in `R`) but I don't have the original (raw) data. I have only the correlation matrix of all the variables. I have seen some questions here about this, but my question continue unsolved. A user gave a parcial solution (http://www.stat.wmich.edu/wang/561/egs/Rcancor.html), but I need the canonical loadings, the percentual of variance in set Y that was explained from set X, and the variables significance. Does anyone here could help me? P.S.: I am a new `R` user. I have experience only on Eviews, GRETL and SPSS (also a little bit in Stata).


title: What happens to adjusted R squared as sample size increases?

text: What effect does sample size have on adjusted R squared values?


title: Calculate odds ratio confidence intervals from plink output?

text: I have output from plink haplotype analysis, however I do not have the raw data. Here, is the output for Haplotype-based association tests with GLMs:               SNP1    SNP2    HAPLOTYPE   F   OR  STAT    P     rs1 rs2 22  0.00992 4.23    61.5    4.43E-15     rs1 rs2 12  0.038   1.02    0.217   0.642     rs1 rs2 21  0.00015 5.22E-10    453 1.77E-100     rs1 rs2 11  0.952   0.762   22.9    1.73E-06      Here is the explanation for each column from plink:                       SNP1    SNP ID of left-most (5') SNP             SNP2    SNP ID of left-most (3') SNP        HAPLOTYPE    Haplotype                 F    Frequency in sample               OR    Estimated odds ratio             STAT    Test statistic (T from Wald test)                P    Asymptotic p-value      Question: based on above output, is it possible to calculate `OR` 95% confidence intervals?


title: Implementing Latent Dirichlet Allocation - notation confusion

text: I am trying to implement LDA using the collapsed Gibbs sampler from http://www.uoguelph.ca/~wdarling/research/papers/TM.pdf the main algorithm is shown below ![enter image description here](http://i.stack.imgur.com/X9OwX.png) I'm a bit confused about the notation in the inner-most loop. n_dk refers to the count of the number of words assigned to topic k in document d, however I'm not sure which document d this is referring to. Is it the document that _word_ (from the next outer loop) is in? Furthermore, the paper does not show how to get the hyperparameters alpha and beta. Should these be guessed and then tuned? Furthermore, I don't understand what the _W_ refers to in the inner-most loop (or the beta without the subscript). Could anyone enlighten me?


title: Does full subset selection regression model building suffer from the same handicaps as stepwise regression?

text: Let's assume $p$ potential predictor variables $X_1,...,X_p$ and a single dependent variable $Y$. Now I evaluate the performance of all possible linear models considering all possible combinations of predictor variables ($2^p-1$). The performance measure(s) could be pretty much any statistic but first comes to mind $R^2$, $F$-statistic and MSE. Based on them I select the "best" model _or_ the top selection which I can check out more closely. Intuitively I would assume this is a great idea--but I read around a bit and came across the infamous concept of "stepwise regression" and how it is considered useless by a lot of ( though apparently not all) statistically trained people. The reason seems to be that the assumed distribution underlying the involved statistics does not hold for the scenario. But stepwise regression is usually described as a slightly different algorithm where you start with a model, adding and removing variables from the model based on a criterion for a statistic. So my question is whether the approach I describe would also be a type of stepwise regression and hence be handicapped by design. On the latter part (if it is SwReg) I would be interested in clarification on where the handicap comes into play and whether it is possible to amend it.


title: Can a -2 Log likelihood be calculated with only one model?

text: I am using the `glmfit` function in MATLAB. The function only returns the deviance and not the log likelihood. I understand that the deviance is basically twice the difference between the log likelihoods of the models but what I don't get is I am only using `glmfit` to create one model, but somehow I am getting a deviance.   * Doesn't calculation of the -2 Log likelihood require 2 models?   * How can the deviance be analysed when there is only one model? Another question I am having is say I did have two models and that I was comparing them using log likelihood test. The null hypothesis would be the first model and the alternative hypothesis would be the second model. After getting the log likelihood test statistic would I check it against the chi squared cdf to determine the p-value? Am I right that if it is less than the alpha level I would reject the null and if it is greater I would fail to reject the null?


title: Probability of an event occurring between two events?

text: 1. I have a data set which looks like this:! http://i.imgur.com/vrDpuQ5.png Histogram: ![http://i.imgur.com/2U5QtQM.png](http://i.stack.imgur.com/eUomJ.png) There are about 10000 entries in total the Outcomes(EventId#) can be either $1,2,3,\ldots,18,19$ so from 1-19.   2. Although it is highly debatable if the events are independent, let's assume for simplicity's sake that they are for now. The reality is that for the smaller EventId's (1,2,3) are more likely to occur after each other and especially after EventId 19. Then again, let's ignore this for now.   3. I am looking specifically for what methodology/formula/principles would allow me to calculate the probability of EventID/Number $n$ to occur between two numbers $x_1$ and $x_2$. Or as can be seen in the first picture of the data set, the probability of the EventId/Number "17" occuring betwen two number/EventId "18". The amount of numbers between $x_1$ and $x_2$ and whether they are located before or after the number $n$, is irrelevant to us, as long as they do not exceed the numbers $n$ and $x_1$ and $x_2$.  So throughout the entire data set we are looking for the probability of events / subseries(?) like these:               2,15,3,6,7,8,2,17,16,2,13,6,7,2,3,1,6,15,         **18,3,6,1,2,9,10,17,14,1,6,16,4,14,18**,2,3,5,1,2      or this:               5,8,2,3,11,15,1,**18,2,3,17,3,11,6,14,16,2,1,6,18**,      2,1,1,9,8,1,13,11,2,7,9,1,2,3,7,2,1,7      So how would I approach solving this problem. What is the solution I am looking for? I have tried looking at almost everything the internet would offer me as a solution, cumulative distribution functions, empirical distribution functions, multivariate distribution, exponential distribution, etc., but from my understanding is that most of these things really just deal with a number/EventID occurring that is between number $x_1$ and Number $x_2$. Which is fine but it doesn't necessarily deal with for example: EventId 18 happening, and then waiting until EventId 17 happens and then waiting again until we have EventId 18. The same obviously could apply to EventId 5 happening then waiting for EventId 10 and then waiting for EventId 15 all the while not allowing any numbers/EventId in between to be higher than 10 and 15. Now assuming I can have / expand this data set into a much much larger data set and I have much more data available, sooner or later the obvious logical conclusion is in fact that EventId 1 2 and 3 tend to follow EventId 19 and 18 more often than for the other EventId's. But I would have to do more analysis on this using Markov chains or some other tool I assume which is why I said for the moment being, let's just assume that the EventId's are independent and random. So yeah, I'm pretty much stuck, if anyone has an idea how to approach this problem then I'd appreciate it, quite a lot.


title: Power of a hypothesis test

text: My question may be pedantic, but I find the syntax used to describe the "power" of a hypothesis test really annoying, and I just want to either correct my understanding of the concept, or understand why it is, seemingly, mislabeled. My understanding of "power" is the probability that a test will reject $H_0$, when $H_a$ is true. Just from this wording, the power of a test seems to be the power _in favor_ of $H_a$, not the power _against_ it. In other words, if $H_a$ is true, you want your test to reject $H_a$. However, my text, _Introduction to the Practice of Statistics_ , used as a text book at our community college, states that this probability is "the power _against_ $H_a$". Why do they call the probability that $H_0$ is rejected when $H_a$ is true the power against $H_0$ instead of the power in favor of $H_a$, or simply the power of $H_a$?


title: Probabilty of getting an open-ended straight draw in 5 cards

text: I play Hold'Em on an internet site and it seems that far too often I see 4 sequential cards being flopped in the full 5 card flop; which is an easy setup for at least one of the players to make a straight. For instance, "3-4-5-6" or "8-9-10-J". I want to know what is probability of getting 4 such sequential cards in 5 cards. It would also be nice to know the same probability as above but without one of cards being an Ace. That's what creates an open-ended straight draw.


title: Can I use glm algorithms to do a multinomial logistic regression?

text: I am using spotfire (S++) for statistical analysis in my project and I have to run multinomial logistic regression for a large data set. I know best algorithm would have been mlogit, but unfortunately that is not available in s++. However, I have an option of using glm algorithm for this regression. I want to clarify two things here: 1.Is my understanding correct that glm also can be used to run Multinomial Logistic Regression?   1. If answer to previous question is yes, then what parameters should be used in glm algo? Thanks,


title: Comparing predictions from different regression equations

text: I have some regression equations, take from various published studies, that predict an individual's mass from its length. Typically only an equation's parameter values, R-squared, standard error of the estimate, and sample size are reported in a manuscript. I'm then using those equations to predict the masses of individuals from a new set of length observations. I'm mainly interested in total mass per unit area, so I'm summing the masses of all individuals and dividing by the area they occupy. Here's a brief example using R code:               # regression equation 1 information     eq1 <- list(b0= 0.9, b1= 3.2, n= 10, r2= 0.984, see= 1.28)          # regression equation 2 information     eq2 <- list(b0= 1.1, b1= 2.8, n= 16, r2= 0.971, see= 1.65)          # new observations     length <- rgamma(100, 4)     area <- 1000          # equation 1 prediction     mass.eq1 <- eq1$b0 + eq1$b1 * length     massPerArea.eq1 <- sum(mass.eq1) / area          # equation 2 prediction     mass.eq2 <- eq2$b0 + eq2$b1 * length     massPerArea.eq2 <- sum(mass.eq2) / area          # compare the two massPerArea predictions...?      Each regression equation of course results in a different final estimate, but how can I determine the uncertainty in those estimates and to what degree they differ statistically? Since I'm summing up individuals, is propagation of error part of that uncertainty? If a direct computation of that uncertainty is possible, that'd would be great, but if there's a solution that requires numerical simulations, that'd be fine also (and actually, I've tried some simulating, but I thought I'd ask first if a direct approach exists). Thanks!


title: How do I test whether differences in populations are significant?

text: I'm working on a large database. I want to explore whether certain batch jobs run in more or less time after a change was than they did before the change. There are about 90 distinct batch jobs. They run according to different timetables, so some of them are called only once a week, and some are called hundred of times per day. I have data for a three-month period before the change was made, and a one-week period since then. Each distinct batch job has run a different number of times before and after the change being investigated. Currently I am using a two-sample z-test, as described here: http://www.stat.ucla.edu/~cochran/stat10/winter/lectures/lect21.html For each distinct batch job, I have sample size, mean run time, and standard deviation, both before and after the change. So I'm calculating a z-value, based on a null hypothesis that the two populations have the same mean and s.d. The resulting z-values range from -35 to +10. My question is, is this the right statistical test for this situation? I believe that the z-values tell me how many standard deviations away from a mean result, given the null hypothesis, my new (post-change) population is. I believe I can use that information to identify the batch jobs that have been most affected by the change (accepting that correlation ! => causation). Am I doing the right test? Have I understood the meaning of the data correctly, and what I can infer from it?


title: Capturing the relationships among model variables based on ANOVA table

text: While studying ANOVA analysis, I meet an example related to the ANOVA table given in the attached figure ![enter image description here](http://i.stack.imgur.com/VuoaG.jpg) Here, > $Y$ stands for change in hemoglobin (%) > > $X_1$ stands for duration of the operation (min) > > $X_2$ stands for blood loss (ml) Based on this ANOVA table, several arguments are made:   * $X_2$ has a significant linear association with $Y$ with or without including $X_1$.   * $X_1$ has a significant linear association with $Y$ after adjusting for linear effects of blood loss X2 on both $X_1$ and $Y$.   * Without adjusting for linear effects of $X_2$, the linear relationship between $Y$ and $X_1$ was not quite significant I am quite confusing on how to generate these three arguments based on the information of ANOVA table? I think it should have connection with those P-values, but do not know how?


title: Testing whether the ratio of two sums is equal to a certain value

text: I have the following hypothesis: $h_0 : \frac {\sum {a_{i}}} {\sum {b_{i}}} = x \\\ h_1 : \frac {\sum {a_{i}}} {\sum {b_{i}}} \neq x$ Where $a_i \propto X b_i$ and $X$ is a random variable. How would I test this hypothesis? Edit: Switched a and b around in the 3rd equation


title: Verifying the assumptions of the paired t-test

text: I have data on the number of cases that customers open (with an engineer) at a computer company each month, across the various technologies that are supported. The assumption is the documents are located on the company's support website / the customers actually view the documents before opening a case. These data include cases that were opened before and after any documentation was created. I am interested in studying the overall effect that the documentation has on case deflection. I am considering using the paired t-test but I have a few questions about verifying the assumptions, and about other reasonable tests to use if the assumptions aren't satisfied. For this example say there are 35 different technologies, and there are 200 supporting documents divided up unequally among the 35 technologies. Below is a an example to make my problem more clear, followed by my questions. **Example:** Say there is a year's worth of data for technology $x_1,~x_2,~...,~x_{35}$ pre / post documentation creation. For instance, technology $x_1$ may have had 568 cases opened prior to documentation creation, then say we noticed that 276 cases were opened for $x_1$ after 46 documents were created for $x_1$. For $x_2$ there were 438 cases that were created prior to documentation creation, and 155 cases opened after 27 documents were created for $x_2$, etc. I want to study the overall effect that the documents as a whole had on case deflection. I listed my questions as bullet points:   * Since one document is not being applied to all of the technologies can I still use the paired t-test? Does it make sense to use the paired t-test in this example?   * One assumption of the paired t-test is that the paired differences are normally distributed. What efficient / easy to use tests that I can utilize to check for normality? I can utilize both R and Excel for the check.   * Another assumption for the paired t-test is that the pairs are independent. How does independence apply to my situation?   * If the paired differences are not normally distributed, then what other test can I use for my problem? I came across the "Wilcoxon signed-rank test", would this be an appropriate test to use? If so, what assumptions must be met to use it / would it apply in my problem?


title: What does "fiducial" mean (in the context of statistics)?

text: When I Google for               "fisher" "fiducial"      ...I sure get a lot of hits, but all the ones I've followed are utterly beyond my comprehension. All these hits do seem to have one thing in common: they are all written for dyed-in-the-wool statisticians, people thoroughly steeped in the theory, practice, history, and lore of statistics. (Hence, none of these accounts bothers to explain or even illustrate what Fisher meant by "fiducial" without resorting to oceans of jargon and/or passing the buck to some classic or other of the mathematical statistics literature.) Well, I don't belong to the select intended audience that could benefit for what I've found on the subject, and maybe this explains why every one of my attempts to understand what Fisher meant by "fiducial" has crashed against a wall of incomprehensible gibberish. Does anyone know of an attempt to explain _to someone who is not a professional statistician_ what Fisher meant by "fiducial"? P.S. I realize that Fisher was a bit of a moving target when it came to pinning down what he meant by "fiducial", but I figure the term must have some "constant core" of meaning, otherwise it could not function (as it clearly does) as terminology that is generally understood within the field.


title: CHAID vs CRT (or CART)

text: I am running a decision tree classification using SPSS on a data set with around 20 predictors (categorical with few categories). CHAID (Chi-squared Automatic Interaction Detection) and CRT/CART (Classification And Regression Trees) are giving me different trees. Can anyone explain the relative merits of CHAID vs CRT? What are the implications of using one method over the other?


title: Statistics wording question

text: I am reading a statistical procedure trying to figure out and understand what's going on. The statement says "Compute the posterior on $\mu$." Does this mean compute $p(\mu)$? Does this mean "Compute the posterior probability of $\mu$." (see how "of" replaces "on")? Immediately before in the reading, $\mu$ is allowed to take on multiple values. Maybe there is a statistical wording/vocabulary convention I'm unaware of that would make this make more sense to me? This is not critical, but I am just trying to understand better. If anyone could help me or point me to a resource that would be great.


title: Sampling distribution of regression coefficients

text: I previously learned about sampling distributions that gave results which were for the estimator, in terms of the unknown parameter. For example, for the sampling distributions of $\hat\beta_0$ and $\hat\beta_1$ in the linear regression model $Y_i = \beta_o + \beta_1 X_i + \varepsilon_i$ $$ \hat{\beta}_0 \sim \mathcal N \left(\beta_0,~\sigma^2\left(\frac{1}{n}+\frac{\bar{x}^2}{S_{xx}}\right)\right) $$ and $$ \hat{\beta}_1 \sim \mathcal N \left(\beta_1,~\frac{\sigma^2}{S_{xx}}\right) $$ where $S_{xx} = \sum_{i=1}^n (x_i^2) -n \bar{x}^2$ But now I have seen the following in a book: > Suppose we fit the model by least squares in the usual way. Consider the > Bayesian posterior distribution, and choose priors so that this is > equivalent to the usual frequentist sampling distribution, that is...... > > $$ \left( \begin{matrix} \beta_0 \\\ \beta_1 \end{matrix} \right) \sim > \mathcal N_2\left[\left(\begin{matrix} \hat{\beta}_1 \\\ \hat{\beta}_2 > \end{matrix} \right),~\hat{\sigma}^2 \left(\begin{matrix} n & > \sum_{i=1}^{n}x_i \\\ \sum_{i=1}^{n}x_i & \sum_{i=1}^{n}x_i^2 \end{matrix} > \right) ^{-1}\right] $$ This is confusing me because:   1. Why do the estimates appear on the left hand side (lhs) of the first 2 expressions, and the right hand side (rhs) of the last expression?   2. Why do the beta hats in the last expression have 1 and 2 subscripts instead of 0 and 1?   3. Are these just different representations of the same thing? If they are, could someone show my how they are equivalent? If not, could someone explain the difference?   4. Is it the case that the last expression is the "inversion" of the first two? Is that why the 2x2 matrix in the last expression is inverted and estimates/parameters are switched from rhs$\leftrightarrow$lhs? If so could someone show me how to get from one to the others?


title: What statistical test should I use?

text: I have 10 different biological features (for example: normal, disease 1 , disease 2, ..., disease 9). Tissue samples from various patients were imaged and 60 imaging parameters were measured. A clinician identified different regions in the imaged tissue to be normal or diseased (1-9). I want to know if the measured parameters can distinguish between the different tissue conditions. I do not know if these parameters are independent of one another. There are thousands of values (for each parameter) assigned to a tissue condition. But this data was obtained by imaging a few tissue samples. What sort of statistical test should I use to determine if the 60 parameters are statistically significantly different between the 10 different features? Is there anything in MATLAB? I was thinking of `anovan.m`, but it's not clear how I should go about using it.


title: Best data series to teach time series econometrics?

text: I will be teaching a short course on time-series econometrics to third-year undergraduates. It will be part theoretical, part applied (students using econometrics software). I am looking for great examples of data series which will fascinate (inspire?!) the students and engage them with time series econometrics. Topics of the course will include:   * Autoregressive models   * Autoregressive distributed lag models   * Non-stationarity (basic unit root testing)   * Engle-Granger type cointegration analysis   * Structural breaks Have you come across great example data sets with interesting properties and interpretations that you can recommend for some of the above applications? I would much appreciate it! Thank you!


title: Standard score (z-transformation) normalization or 0-1 normalization, which one is better for k-NN?

text: I am not much experienced in data mining, but I know that I should normalize my data before running k-NN classifier on it, to have reasonable results. But I found out, that there are many methods to do that, and would like to know, which one would you suggest? Which one is better for k-NN?               standard score (z-transformation) normalization       or               0-1 normalization       example - result of z-transformation:                statistics              range     avg = 0 +/- 1   [-0.651 ; 7.158]     avg = 0 +/- 1   [-1.156 ; 3.197]     avg = 0 +/- 1   [-0.023 ; 74.710]     avg = 0 +/- 1   [-0.141 ; 21.766]     avg = 0 +/- 1   [-0.017 ; 75.289]     avg = 0 +/- 1   [-0.094 ; 49.605]     avg = 0 +/- 1   [-0.057 ; 42.350]     avg = 0 +/- 1   [-0.065 ; 57.510]     avg = 0 +/- 1   [-0.054 ; 50.810]     avg = 0 +/- 1   [-0.056 ; 65.837]     avg = 0 +/- 1   [-0.065 ; 40.408]     avg = 0 +/- 1   [-0.150 ; 11.658]     avg = 0 +/- 1   [-0.442 ; 6.869]     avg = 0 +/- 1   [-0.196 ; 52.817]     avg = 0 +/- 1   [-0.349 ; 27.916]     avg = 0 +/- 0   [0.000 ; 0.000]     avg = 0 +/- 1   [-2.065 ; 16.025]     avg = 0 +/- 1   [-0.356 ; 42.274]     avg = 0 +/- 1   [-3.532 ; 1.142]               statistics              range     avg = 0.083 +/- 0.128   [0.000 ; 1.000]     avg = 0.266 +/- 0.230   [0.000 ; 1.000]     avg = 0.000 +/- 0.013   [0.000 ; 1.000]     avg = 0.006 +/- 0.046   [0.000 ; 1.000]     avg = 0.000 +/- 0.013   [0.000 ; 1.000]     avg = 0.002 +/- 0.020   [0.000 ; 1.000]     avg = 0.001 +/- 0.024   [0.000 ; 1.000]     avg = 0.001 +/- 0.017   [0.000 ; 1.000]     avg = 0.001 +/- 0.020   [0.000 ; 1.000]     avg = 0.001 +/- 0.015   [0.000 ; 1.000]     avg = 0.002 +/- 0.025   [0.000 ; 1.000]     avg = 0.013 +/- 0.085   [0.000 ; 1.000]     avg = 0.061 +/- 0.137   [0.000 ; 1.000]     avg = 0.004 +/- 0.019   [0.000 ; 1.000]     avg = 0.012 +/- 0.035   [0.000 ; 1.000]     avg = 0.010 +/- 0.028   [0.000 ; 1.000]     avg = 0.114 +/- 0.055   [0.000 ; 1.000]     avg = 0.008 +/- 0.023   [0.000 ; 1.000]     avg = 0.756 +/- 0.214   [0.000 ; 1.000]


title: How do we interpret the results of regression for log transformed data containing negative values, $\log(Y-\min(Y)+1)$?

text: I have a question regarding interpretation of regression results based on data where I have some negative values. Since the residuals were positively skewed I needed to log-transform my data, and because of negative values my log transformation looked like $\ln(Y-\min(Y)+1)$. Normally, if we have $\ln$ in both sides of the equation, then our $b$ coefficient means a percentage change of DV for 1 % change of IV. But I have big doubts whether the same is valid for such log-transformed data, where I added a constant to my DV in order to make the minimum value equal to 1. I can't also find any literature on this topic. Does anyone know for sure what would be the correct interpretation in such a situation?


title: How can I calculate the R-square value of an MDS to assess the fit of the model to the data in R?

text: I would like to do Multidimensional Scaling (MDS) using `cmdscale()` in R. I have read that it is useful to try out how many dimensions are suitable for the data by trying different values of _k_ , and then seeing what proportion of variance is accounted for in the MDS result by looking at the R-square value. R-square values smaller than 0.6 are generally found to be acceptable for a good fit between the data and the number of dimensions. However, how do I calculate R-square from an MDS generated by `cmdscale()`?


title: Logistic regression with categorical predictors, do log-odds differ from 0?

text: (My title is probably not very clear, but I don't know how to state it more clearly. So feel free to edit it. Thanks!) I have a 3 by 2 design, with a total of 6 conditions. The outcome is binary (0 or 1). Below is a sample dataset generated in `R`:               set.seed(2)     mockdata<-data.frame(outcome=sample(1:0, 48, prob=c(0.5, 0.5), replace=TRUE),                          f1=rep(letters[1:2], each=24),                           f2=rep(letters[1:3], each=8))          head(mockdata)     #  outcome f1 f2     #1       0  a  a     #2       1  a  a     #3       1  a  a     #4       0  a  a     #5       1  a  a     #6       1  a  a      _One of the things I would like to look at is whether the log-odds of the outcome for each of the 6 conditions is significantly different from 0._ I can create a new condition variable as follows:               mockdata$f12 <- paste(mockdata$f1, mockdata$f2, sep=".")      then, I can do logistic regression using the newly created variable (see below for output). The intercept below tells me that for the condition that is treated as the baseline condition, the log-odds is not significant different 0. **My questions are:** **(1). to check the other conditions, should I simply change the baseline condition, and after testing all 6 conditions, I adjust the p-values accordingly?** **(2). Are there better ways of testing what I want to test?**               summary(glm(outcome ~f12, family="binomial", data=mockdata))          Call:     glm(formula = outcome ~ f12, family = "binomial", data = mockdata)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -1.6651  -1.1774  -0.5168   1.0215   2.0393            Coefficients:                 Estimate Std. Error z value Pr(>|z|)       (Intercept)  -0.1412     0.3281  -0.430   0.6669       f121          0.6520     0.6806   0.958   0.3380       f122          0.1412     0.6641   0.213   0.8316       f123         -0.3696     0.6806  -0.543   0.5871       f124         -1.8047     0.9325  -1.935   0.0530 .     f125          1.2398     0.7430   1.669   0.0952 .     ---      **EDIT** I would also like to check how the log-odds differ amongst the different conditions. For that, I was gonna just run a regular logistic regression with the main effects of f1 and f2, and the interaction of the two, and conduct additional multiple comparisons dependent on the kind of result I get from the omnibus test.