title: Effects of Re-coding a regression

text: I have been analysing a number of likert scale variables using spatial auto- regressive models. However, my question I think can be applied to an OLS model. After analyzing the data I obtained some odd coefficients, which led me to realize that some of the variables were coded inversely. That is, 1 was very important and the highest number was not very important. So I re-coded the response so that the larger the number the more important it is. I had assumed that this would only change the coefficients direction, and maybe slightly effect the magnitude, however I have received some dramatically different results. My question is then: Should I have expected to see such a dramatic change in the results? Or am I correct in being suspicious of my re-coding?


title: Outputting only p-value for ANOVA in R

text: I'm doing some basic one-way ANOVA in R. It would be useful to be able to output only the p-value associated with a _specific_ term without producing the whole ANOVA table (unrelated reason). I need the value as a basic numeric and not an ANOVA object For example, if I'm looking at `anova(lm(yield~variety+block))` and just want the p-value for "variety", is there a straightforward way? Thanks!


title: Is there a population size limit for using the R package Demography when doing stochastic population forecasting

text: Is it possible to do stochastic population forecasting with the R package Demography when the population size is between 10k and 20k? I have data from year 1985 to 2012 for population, deaths and births. What do you do when almost half of the death rates are 0? From ages 1 to 30 nearly no one dies in a so small population. Is there a method for imputing some plausible values for the zero's in the death and birth rates? Is it possible to, for example, multiply all population numbers with 1000 and impute some plausible values for rates (maybe borrow these from a similar population), and then divide all the simulated numbers in output with 1000? I have tried to do some forecasting with the R package demography, but I get so many warnings and weird results, I can't get it right. Is there a population size limit for using the `demography` package? Do anyone have any suggestions for how to solve the problem of stochastic forecasting with small populations?


title: How to analyze repeated measurements in two matched but unequal samples

text: For a manuscript intended for submission to a medical (or psychological) journal, I have cognitive data from 16 patients and 32 age- and gender-matched healthy controls at four time points (pre-treatment T0, and three post- treatment time points T1, T2, and T3 in patients; in controls at the same intervals), i.e., I have two healthy controls matched with each patient. The cognitive measure, though technically discrete, can for most purposes be thought of as normally distributed. Example of what the data looks like:               subject gender age  T0  T1  T2  T3     p1      male   12   30  39  44  45     p2      female 15   37  38  40  36     ...     ...    ...  ... ... ... ...     ca1     male   12   38  41  52  48     ca2     female 15   43  42  40  41     ...     ...    ...  ... ... ... ...     cb1     male   12   29  39  49  61     cb2     female 15   49  52  53  59     ...     ...    ...  ... ... ... ...      I want to know: 1\. Do patients differ from controls before treatment (T0)? 2a. Do patients differ from controls after treatment? 2b. Do pre-existing differences as in 1, if any, change? However, I know that the cognitive measure is (normally, i.e., in healthy controls) significantly predicted by   a. age, gender, and their interaction   b. repeated measurement I would therefore like to retain the matching in the analysis. However, I am not sure how to achieve this. Would it be feasible to perform a repeated measures analysis of variance with two "within-subject" factors (time: T0,T1,T2,T3, and subject type: patient, controla, controlb) and then use contrasts to achieve the desired grouping? If so, how do I set this up properly in R? Another option I came up with is to calculate a set of 32 difference scores for each assessment (p1T0-ca1T0,...,p16T0-ca16T0,p1T0-cb1T0,...,p16T0-cb16T0) and then use one-sample t-tests to test for the significance of the average difference (although that may require explicit correction for multiple comparisons). However, I am not certain that is statistically valid (at least in terms of the degrees of freedom applied, normally 32-1 = 31). Is the anova technique feasible and proper? Should I go with the simple t-test method, and in that case, would it be correct to use 31 degrees of freedom? Or do I need another class of methods entirely?


title: Estimate quantile of value in a vector

text: I have a set of real numbers. I need to estimate the quantile of a new number. Is there any clean way to do this in R? in general? I hope this is not ultra-trivial ;-) Much appreciated for your response. PK


title: Dummy variables and likert scale

text: So I have two related questions. I am performing a logistic regression model and want to know how to enter these variables. 1) For one survey I have questions and the responses are Yes, Don't know, and No. Can I just create a dummy variable for Yes (1) with Don't know and No coded as 0 or do I need to exclude the Don't know's (this significantly cuts down my sample size). 2) I have a likert scale for another survey (1 strongly agree, 3 neither, 5 strongly disagree) what is the best way to enter this into the model? Just as categorical data. Some people have argued to use as continuous or again to create dummy variables.


title: Time series and anomaly detection

text: I would like to setup up an algorithm for detecting an anomaly in time series, and I plan to use clustering for that.   * Why should I use a distance matrix for clustering and not the raw time series data?,    * For the detection of the anomaly, I will use density-based clustering, an algorithm as DBscan, so would that work in this case? Is there an online version for streaming data?   * I would like to detect the anomaly before it happens, so , would using a trend detection algorithm (ARIMA) be a good choice?


title: Measures of model complexity

text: How can we compare complexity of two models with the same number of parameters? **Edit 09/19** : To clarify, model complexity is a measure of how hard it is to learn from limited data. When two models fit existing data equally well, a model with lower complexity will give lower error on future data. When approximations are used, this may technically not always be true, but that's OK if it tends to be true in practice. Various approximations give different complexity measures


title: What is the most appropriate way to combine three NTP server results?

text: I need to query an internet time server, in order to assess the staleness of streamed data I am receiving on another channel. I periodically make 3 successive queries, then use the combined error to adjust my system clock until the next set of 3 queries. Each query returns an accurate time from the remote server. I adjust this time according to the time taken to get the reply, assuming symmetric network propagation delay. Typical results might look like this: 21ms, 18ms, 23ms 28ms, 31ms, 76ms At the moment I'm using a mean average, but I'm not convinced this is the best way. In the first case, the mean average seems reasonable. But in the second case, it seems the third figure was affected by a transient asymmetric network delay, and should probably be disregarded. How can I use stats to help me here?


title: Creating tree splits in Hough forest

text: I am working on my thesis and working machine learning. The overall problem is Detecting and Recognition of Road Inventory. For my research part I am looking into decision trees, especially Hough forests. There are some areas where I need more information and if anyone else have worked with this I might save some time and get more time to do new scientific research. So here are a few questions related to the area, that I would like to discuss. (I will reference this question in my thesis if anyone bring up some scientific research). Based on Hough Forests for object detection, tracking and action recognition by Gall, Juergen et al. and other related articles using random forests (forests of decisions trees) most of what I have come across have described and used training splits based on checking if the difference of two features (etc. the x-gradient of the p and q pixel position) in a image patch are lower then a split value tau. Then the best split are found where the information gain is largest, (i.e. splitting the data into the two most distinctive subgroups). A few years earlier then these articles, from what I get it all origin from the work of Breiman. The implementation is available in OpenCV.   1. Have there been any work in the two methods, the OpenCV implementation and the outlined idea in Gall et al. The overall difference is the way the split are done. Gall et al compare two features at p and q, OpenCV select one feature and find the best split for it. Is one of the methods to prefer?   2. The hole scale part of Gall et al. work. From what I get, all positive objects are scaled such longest spatial dimension to a unit size around 100 pixel. I get the reason for doing it, but when taking into fact that the test image are not scaled at all, will this not mean they are trying to match an image patch of size X with dictionary image patches of another size. Cant this be done in a better way? Might create follow up questions as my work continues, but I would like to discuss these thins, as I find them interesting for my project and I am sure there are people out there who worked tons of hours in similar topics.


title: Estimating parameters of a normal distribution from noisy observation of samples

text: Suppose I have a number of samples drawn from a normal distribution $x_i \sim \mathcal{N}(\mu,C)$ with $i = 1 \dots n$. I can make observations $z_i = x_i + e_i$ for those samples which are perturbed by noise $e_i$ of known characteristic $e_i \sim \mathcal{N}(0,\Sigma_i)$. Note that the distribution of the observation noise is different for each $z_i$. I would like to estimate $\mu$ and $C$ given the $z_i$. Since $e_i$ is unbiased $\mu$ should just be $\mathbb{E}[z_i]$. In case $\Sigma_i = D$ for all $i$, $C = \operatorname{Var}(z_i) - D$ assuming independence between the noise. Would it be safe to say that $C = \operatorname{Var}(z_i) - \mathbb{E}[\Sigma_i]$ in the general case?


title: Checking ANOVA assumptions

text: A few months ago I posted a question about homoscedascity tests in R on SO, and Ian Fellows answered that (I'll paraphrase his answer very loosely) "homoscedascity tests are not a good tool when testing the goodness of fit of your model. With small samples, you don't have enough power to detect departures from homoscedascity, while with big samples you have "plenty of power", so you are more likely to screen even trivial departures from equality". His great answer came as a slap in my face. I used to check normality and homoscedascity assumptions each time I ran ANOVA. Now you lads hop in: What is, in your opinion, best practice when checking ANOVA assumptions?


title: Estimating parameters for probit multiplied by something

text: Assume you have a model of the form $y = x_1\Phi(\beta_0 + \beta_2x_2 +\ldots+ \beta_nx_n)+u$ where $y_i\in(0,x_{1i}]$ and $\Phi$ is the probit function. How can we estimate $\beta$s in e.g. Stata?


title: Variations in time series

text: I'm trying to do some really simple Time Series work. I started off looking at monthly counts but the more I thought about it (and the more I read posts by IrishStat) the more I realised my view of this was fundamentally flawed. My main reason for thinking this is the plethora of timing issues there are with using monthly aggregate figures:   * Varying month length   * Varying numbers of weekdays   * Weekends   * Holidays But what gets me is that this must be a really common issue, I mean practically every big company out there must have an executive report that includes monthly figures. How does one balance for the discrepancies between the months you are comparing? I'm using R to slowly crawl through this stuff and currently I'm just looking at cross correlations with lag. I've looked at decomposition and how to adjust for seasonality but how does one adjust for dates that you know will have an effect on business activity? Is there a way to integrate a list of dates that should be treated differently or somehow compensate for the number of weekend/holiday days in a month (or in a week, or in a year)? A general explanation of the concept would be much appreciated. **Edited to make the question more succinct:** I am looking at a time series of monthly counts. 2 examples from this time series would be:   * February 2012 is 29 days long, has 8 weekend days and zero public holidays   * May 2013 is 31 days long, also has 8 weekend days but also has 2 public holidays How would I compare them without having my conclusions adversely affected by their differences?   * I could divide by days in month but this doesn't take into account holidays   * I could divide by working days in month, but this doesn't differentiate between a holiday and a weekend.


title: What could be the described method for treating time-outs in algorithm benchmarks

text: As a computer scientist, you often face the problem to empirically analyse the improvements in run time of some algorithm. I stumbled over the following text in a paper: > These run times are averages over 10 runs per benchmark, and account for > time-outs using maximum likelihood estimation (MLE) [cite]. With 90% > confidence level, 71% of those averages are accurate to within 25%. The cite is simply a book on statistics "Mathematical Statistics and Data Analysis, Rice" without page numbers :( ## To explain the setting Suppose you have a set of algorithms $A$. Every algorithm is applicable to every benchmark problem in your problem set $P$. This gives $|A||P|$ possible configurations. For randomized algorithms, you repeat each configuration $n$ times. This makes a total of $n|A||P|$ experiments. For each experiment, you measure the time it takes the algorithm to calculate the solution. But for practical reasons, you set a time-out $t_{max}$. So you end up with a set of run-times with na values meaning "larger than $t_{max}$". The cite is from a caption plotting number of solved instances over time. So I assume they calculated some kind of "average value" for each algorithm/problem combination (over 10 different random seeds), accounting for the time-outs.


title: What are good resources/criteria for judging human bias in data collection?

text: I've just been given a stack of polling data to analyse. Some of the questions are obviously leading or present subtle incentives (for the poller or polled) for specific answers. Of other questions I'm not so sure but I have some doubts. I'm also starting to question other factors about how the poll was conducted (environment, privacy, etc.). I'd like to be able to present the analysis with the unreliable data discounted according to some recognised standard. What are actionable rules or methods for assessing whether results are junk based on human factors in data collection? I'd gladly accept specific criteria, or links to resources such as books or web pages on the subject. As I'm not a professional statistician, even being told what terms are used for this field of study and it's concerns would help a lot (though maybe best added as comments).


title: Hypothesis testing for normal distribution

text: This question seems easy but I am not getting it. Let $X_{1},X_{2}, · · · ,X_{15}$ be a random sample from a normal population $X ∼ N(0, σ^2)$. Find a best critical region of size $α = 0.05$ for testing $H_{0} : σ_2 = 1$ versus $H_{2} : σ_{2} = 4$.


title: How many tests do I have to run for this case?

text: I have developed a method to estimate a certain value called "flow" from 10% of the total data received (the total data may change every time) but of course there is an error associated with my method. My data comes from simulation that I can run as many times that I want. with how many 10% samples do I have to test my method to have a reliable estimate of the error of my method.


title: Multiple comparisons with t-test: one group with different variance

text: If I have 3 groups that I want to compare with a t-test (correcting for multiple comparisons afterwards) but one of the groups has different variance than the other two, should I use a Welch correction for all comparisons or just for the comparison involving the group with unequal variance?


title: How can I interpret binary logistic table?

text: I´m beginner with SPSS and I have on problem on interpreting binary logistic results. So I have this table: Variables in the Equation                                         B      S.E.  Wald     df   Sig.   Exp(B)     Step 1ais_first_timer   -1.582  .129  149.593   1   .000     .206         trip_duration       -.301   .047  41.114    1   .000     .740         back_to_country      .583   .130  20.030    1   .000    1.791         back_to_event       1.275   .255  25.059    1   .000    3.578         visiting_places     -.179   .120   2.203    1   .138     .836         Constant            1.072   .145  54.862    1   .000    2.922      So my problem is how can I write an analysis based on that? I know I have to watch Sig and B. But what do they show me? My main problems is what does B shows me? For example if my B is -1.582 and sig is 0.000 and my dependent variable is nationality and covariates are, as you can see in the table, first timer, trip duration etc. Can I say watching this table, for example English people (0) are more likely to be first time visitors than Americans (1), and this is statistically important? Or can i say that English people have longer vacation than Americans Or i can`t say that? I hope somebody can help me :)


title: Compare 2 regression lines in R

text: I have the following data frame               structure(list(Chi = structure(c(1L, 1L, 5L, 5L, 6L, 9L, 9L,      12L, 13L, 14L, 14L, 16L, 16L, 19L, 19L, 20L, 20L, 23L, 24L, 24L,      26L, 26L, 31L, 31L, 33L, 33L, 36L, 37L, 37L, 40L, 40L, 43L, 43L,      44L, 44L, 45L, 46L, 47L, 47L, 48L, 48L, 52L, 52L, 54L, 54L, 55L,      55L, 56L, 59L, 59L, 61L, 61L, 63L, 63L, 64L, 64L, 65L, 65L, 69L,      69L, 70L, 70L, 71L, 71L, 72L, 72L, 75L, 75L, 76L, 76L, 77L, 77L,      79L, 79L, 86L, 86L, 87L, 87L, 88L, 88L, 91L, 91L, 92L, 92L, 93L,      95L, 96L, 96L, 97L, 97L, 98L, 98L, 99L, 99L, 100L, 100L, 101L,      101L, 103L, 103L, 104L, 104L, 107L, 108L, 108L, 112L, 112L, 113L,      116L, 116L, 117L, 120L, 125L, 125L, 127L, 127L, 129L, 129L, 130L,      131L, 131L, 132L, 132L, 134L, 134L, 135L, 135L, 136L, 136L, 139L,      141L, 141L, 143L, 144L, 144L, 145L, 145L, 146L, 150L, 150L, 151L,      151L, 153L, 153L, 155L, 155L, 157L, 162L, 162L, 163L, 163L, 164L,      164L, 167L, 167L, 168L, 169L, 169L, 171L, 171L, 172L, 172L, 174L,      174L, 175L, 175L, 177L, 177L, 180L, 180L, 183L, 187L, 27L, 83L,      83L, 165L, 165L, 85L, 85L, 156L, 156L, 17L, 17L, 123L, 123L,      124L, 124L, 57L, 57L, 42L, 42L, 159L, 159L, 38L, 38L, 82L, 82L,      41L, 41L, 142L, 142L), .Label = c("0106610856", "0107470802",      "0108490513", "0108590534", "0109480651", "0111290260", "0111410339",      "0201390418", "0207570604", "0208360352", "0212323105", "0212380362",      "0301310432", "0302705635", "0303450495", "0304260266", "0304440574",      "0305280546", "0305380338", "0305381393", "0305510576", "0305542214",      "0308610733", "0309370345", "0309665035", "0310380545", "0403320259",      "0403360374", "0404360343", "0406270198", "0501451137", "0504460676",      "0511310366", "0605270511", "0605340560", "0605410461", "0605410585",      "0606260684", "0606270353", "0609360507", "0702520535", "0702570818",      "0705430421", "0710380364", "0801330378", "0801430275", "0802320430",      "0803510802", "0805390383", "0806560533", "0809430460", "0902380354",      "0904340252", "0904370445", "0906340403", "0907380379", "0909415420",      "0910300100", "0911430253", "1001270460", "1001360389", "1002455294",      "1005280487", "1006330445", "1009350447", "1010375156", "1011270447",      "1012350312", "1012400441", "1102570648", "1105450589", "1106230566",      "1106330587", "1204530475", "1206350342", "1208330373", "1209280345",      "1209400502", "1209400561", "1210380536", "1302240455", "1305751256",      "1306370353", "1307260470", "1310340250", "1312430613", "1312440597",      "1312690593", "1404430512", "1404530479", "1405330376", "1406310360",      "1406350419", "1406430439", "1408460602", "1412360366", "1502385236",      "1503370488", "1503470628", "1503660400", "1506390447", "1508340196",      "1510340688", "1510440453", "1603310622", "1604440376", "1606370014",      "1609650549", "1610345304", "1610345304x", "1612300367", "1702330397",      "1704330181", "1706330316", "1712560522", "1802340270", "1804310336",      "1808430417", "1810400244", "1902340299", "1902610679", "1905360355",      "1906320438", "1906390525", "1909310514", "1912460408", "2002440204",      "2004350288", "2007350203", "2009400364", "2009460669", "2011410428",      "2011500524", "2103335236", "2109370262", "2112290355", "2201330484",      "2201600686", "2203290471", "2203406259", "2205430513", "2207340473",      "2208340396", "2303430410", "2303530717", "2308290390", "2309420506",      "2310370398", "2310370398.0", "2312280310", "2404436295", "2406640663",      "2411420404", "2501520858", "2505330239", "2505380376", "2511320428",      "2511320436", "2511360306", "2601490470", "2601520566", "2608450598",      "2611400237", "2701470625", "2702230407", "2702340342", "2703470916",      "2704380538", "2709250586", "2712350545", "2712541146", "2805310438",      "2805350472", "2807360475", "2807480594", "2809325316", "2809470634",      "2902400411", "2903350442", "2905330376", "2906450480", "2910240363",      "3004510529", "3007230195", "3012410333", "3107440299", "3108350420"     ), class = "factor"), Sex = structure(c(2L, 2L, 2L, 2L, 1L, 1L,      1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,      2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,      2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,      1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L,      1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L,      2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,      2L, 2L, 2L), .Label = c("F", "M"), class = "factor"), Age = c(50L,      50L, 63L, 63L, 83L, 55L, 55L, 72L, 81L, 42L, 42L, 86L, 86L, 74L,      74L, 74L, 74L, 50L, 74L, 74L, 73L, 73L, 67L, 67L, 79L, 79L, 71L,      70L, 70L, 75L, 75L, 68L, 68L, 73L, 73L, 79L, 69L, 79L, 79L, 61L,      61L, 74L, 74L, 74L, 74L, 77L, 77L, 73L, 68L, 68L, 76L, 76L, 84L,      84L, 78L, 78L, 77L, 77L, 71L, 71L, 55L, 55L, 67L, 67L, 88L, 88L,      77L, 77L, 78L, 78L, 84L, 84L, 71L, 71L, 69L, 69L, 67L, 67L, 41L,      41L, 78L, 78L, 80L, 80L, 76L, 66L, 76L, 76L, 73L, 73L, 74L, 74L,      64L, 64L, 46L, 46L, 72L, 72L, 78L, 78L, 67L, 67L, 74L, 47L, 47L,      79L, 79L, 79L, 78L, 78L, 81L, 77L, 79L, 79L, 67L, 67L, 76L, 76L,      70L, 64L, 64L, 70L, 70L, 79L, 79L, 74L, 74L, 82L, 82L, 83L, 69L,      69L, 76L, 69L, 69L, 58L, 58L, 83L, 83L, 83L, 68L, 68L, 69L, 69L,      79L, 79L, 79L, 66L, 66L, 70L, 70L, 65L, 65L, 65L, 65L, 72L, 87L,      87L, 57L, 57L, 80L, 80L, 76L, 76L, 63L, 63L, 64L, 64L, 78L, 78L,      60L, 76L, 80L, 75L, 75L, 90L, 90L, 78L, 78L, 74L, 74L, 69L, 69L,      80L, 80L, 73L, 73L, 71L, 71L, 56L, 56L, 76L, 76L, 87L, 87L, 38L,      38L, 61L, 61L, 78L, 78L), SBR = c(12.061, 11.447, 9.403, 9.136,      9.747, 8.648, 7.934, 7.914, 9.349, 11.224, 10.433, 4.897, 5.823,      8.683, 8.692, 13.018, 13.386, 7.817, 7.384, 7.518, 11.091, 11.028,      8.372, 8.497, 10.751, 10.488, 4.347, 2.593, 2.203, 6.461, 7.272,      4.581, 4.593, 10.31, 9.004, 10.362, 10.307, 9.266, 10.163, 9.24,      8.732, 8.449, 7.823, 10.427, 10.669, 8.695, 8.729, 8.653, 12.299,      12.158, 11.748, 11.19, 8.431, 8.717, 8.253, 8.412, 6.911, 6.805,      9.468, 11.413, 6.603, 7.697, 7.762, 7.097, 10.607, 8.162, 5.419,      5.575, 7.007, 6.974, 8.708, 8.419, 9.47, 8.42, 8.229, 8.027,      5.294, 4.628, 11.475, 10.328, 7.905, 8.491, 10.724, 9.02, 9.095,      5.754, 9.805, 7.332, 6.669, 5.118, 12.443, 11.972, 13.309, 13.906,      14.963, 15.119, 6.465, 6.38, 6.949, 6.064, 6.541, 6.648, 3.542,      11.148, 11.918, 9.743, 9.795, 6.103, 6.025, 3.917, 7.304, 7.628,      8.092, 7.347, 9.051, 8.206, 10.697, 10.286, 4.564, 10.62, 9.84,      9.105, 7.998, 6.437, 5.707, 6.949, 6.315, 6.165, 6.68, 8.86,      8.326, 8.6, 7.776, 5.193, 5.456, 11.864, 11.381, 6.385, 10.972,      9.87, 9.645, 7.738, 10.096, 9.667, 9.687, 8.255, 4.606, 8.738,      8.519, 7.002, 6.288, 10.425, 10.303, 8.278, 8.342, 6.657, 6.111,      5.928, 13.06, 12.747, 5.545, 5.845, 9.338, 9.534, 9.635, 8.716,      7.765, 7.254, 7.517, 7.317, 7.335, 5.628, 4.864, 7.1, 7.02, 6.734,      5.622, 7.167, 7.391, 6.443, 6.874, 8.373, 7.573, 5.701, 6.355,      6.884, 6.296, 9.097, 9.645, 7.068, 7.252, 6, 5.794, 8.074, 9.267,      12.584, 10.723, 9.39, 9.165, 9.635, 8.814), Diagnosis = structure(c(2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("A", "N"), class = "factor"),          fit = c(10.1654358296645, 10.1654358296645, 9.07109655193284,          9.07109655193284, 7.38749766311491, 9.74453610746002, 9.74453610746002,          8.31347705196477, 7.5558575519967, 10.8388753851917, 10.8388753851917,          7.13495782979222, 7.13495782979222, 8.14511716308298, 8.14511716308298,          8.14511716308298, 8.14511716308298, 10.1654358296645, 8.14511716308298,          8.14511716308298, 8.22929710752388, 8.22929710752388, 8.73437677416926,          8.73437677416926, 7.7242174408785, 7.7242174408785, 8.39765699640567,          8.48183694084657, 8.48183694084657, 8.06093721864208, 8.06093721864208,          8.65019682972836, 8.65019682972836, 8.22929710752388, 8.22929710752388,          7.7242174408785, 8.56601688528746, 7.7242174408785, 7.7242174408785,          9.23945644081464, 9.23945644081464, 8.14511716308298, 8.14511716308298,          8.14511716308298, 8.14511716308298, 7.89257732976029, 7.89257732976029,          8.22929710752388, 8.65019682972836, 8.65019682972836, 7.97675727420119,          7.97675727420119, 7.30331771867401, 7.30331771867401, 7.80839738531939,          7.80839738531939, 7.89257732976029, 7.89257732976029, 8.39765699640567,          8.39765699640567, 9.74453610746002, 9.74453610746002, 8.73437677416926,          8.73437677416926, 6.96659794091043, 6.96659794091043, 7.89257732976029,          7.89257732976029, 7.80839738531939, 7.80839738531939, 7.30331771867401,          7.30331771867401, 8.39765699640567, 8.39765699640567, 8.56601688528746,          8.56601688528746, 8.73437677416926, 8.73437677416926, 10.9230553296326,          10.9230553296326, 7.80839738531939, 7.80839738531939, 7.6400374964376,          7.6400374964376, 7.97675727420119, 8.81855671861015, 7.97675727420119,          7.97675727420119, 8.22929710752388, 8.22929710752388, 8.14511716308298,          8.14511716308298, 8.98691660749195, 8.98691660749195, 10.5021556074281,          10.5021556074281, 8.31347705196477, 8.31347705196477, 7.80839738531939,          7.80839738531939, 8.73437677416926, 8.73437677416926, 8.14511716308298,          10.4179756629872, 10.4179756629872, 7.7242174408785, 7.7242174408785,          7.7242174408785, 7.80839738531939, 7.80839738531939, 7.5558575519967,          7.89257732976029, 7.7242174408785, 7.7242174408785, 8.73437677416926,          8.73437677416926, 7.97675727420119, 7.97675727420119, 8.48183694084657,          8.98691660749195, 8.98691660749195, 8.48183694084657, 8.48183694084657,          7.7242174408785, 7.7242174408785, 8.14511716308298, 8.14511716308298,          7.47167760755581, 7.47167760755581, 7.38749766311491, 8.56601688528746,          8.56601688528746, 7.97675727420119, 8.56601688528746, 8.56601688528746,          9.49199627413733, 9.49199627413733, 7.38749766311491, 7.38749766311491,          7.38749766311491, 8.65019682972836, 8.65019682972836, 8.56601688528746,          8.56601688528746, 7.7242174408785, 7.7242174408785, 7.7242174408785,          8.81855671861015, 8.81855671861015, 8.48183694084657, 8.48183694084657,          8.90273666305105, 8.90273666305105, 8.90273666305105, 8.90273666305105,          8.31347705196477, 7.05077788535132, 7.05077788535132, 9.57617621857822,          9.57617621857822, 7.6400374964376, 7.6400374964376, 7.97675727420119,          7.97675727420119, 9.07109655193284, 9.07109655193284, 8.98691660749195,          8.98691660749195, 7.80839738531939, 7.80839738531939, 9.32363638525553,          7.97675727420119, 7.6400374964376, 8.06093721864208, 8.06093721864208,          6.79823805202863, 6.79823805202863, 7.80839738531939, 7.80839738531939,          8.14511716308298, 8.14511716308298, 8.56601688528746, 8.56601688528746,          7.6400374964376, 7.6400374964376, 8.22929710752388, 8.22929710752388,          8.39765699640567, 8.39765699640567, 9.66035616301912, 9.66035616301912,          7.97675727420119, 7.97675727420119, 7.05077788535132, 7.05077788535132,          11.1755951629553, 11.1755951629553, 9.23945644081464, 9.23945644081464,          7.80839738531939, 7.80839738531939), lwr = c(5.90999794584117,          5.90999794584117, 4.85411038352648, 4.85411038352648, 3.16383967274129,          5.5078318233643, 5.5078318233643, 4.10340735387365, 3.33646831235876,          6.54339594841324, 6.54339594841324, 2.90340734886211, 2.90340734886211,          3.93437871118446, 3.93437871118446, 3.93437871118446, 3.93437871118446,          5.90999794584117, 3.93437871118446, 3.93437871118446, 4.01899333351728,          4.01899333351728, 4.52246825860982, 4.52246825860982, 3.50829974380329,          3.50829974380329, 4.18762073879345, 4.27163348349591, 4.27163348349591,          3.84956354899072, 3.84956354899072, 4.43905717652239, 4.43905717652239,          4.01899333351728, 4.01899333351728, 3.50829974380329, 4.35544561188301,          3.50829974380329, 3.50829974380329, 5.0187327527967, 5.0187327527967,          3.93437871118446, 3.93437871118446, 3.93437871118446, 3.93437871118446,          3.67933199647321, 3.67933199647321, 4.01899333351728, 4.43905717652239,          4.43905717652239, 3.76454793766266, 3.76454793766266, 3.07722713717148,          3.07722713717148, 3.59391587315661, 3.59391587315661, 3.67933199647321,          3.67933199647321, 4.18762073879345, 4.18762073879345, 5.5078318233643,          5.5078318233643, 4.52246825860982, 4.52246825860982, 2.72879720476565,          2.72879720476565, 3.67933199647321, 3.67933199647321, 3.59391587315661,          3.59391587315661, 3.07722713717148, 3.07722713717148, 4.18762073879345,          4.18762073879345, 4.35544561188301, 4.35544561188301, 4.52246825860982,          4.52246825860982, 6.62171309179384, 6.62171309179384, 3.59391587315661,          3.59391587315661, 3.42248381273475, 3.42248381273475, 3.76454793766266,          4.60567896791123, 3.76454793766266, 3.76454793766266, 4.01899333351728,          4.01899333351728, 3.93437871118446, 3.93437871118446, 4.77149984958079,          4.77149984958079, 6.22823093408572, 6.22823093408572, 4.10340735387365,          4.10340735387365, 3.59391587315661, 3.59391587315661, 4.52246825860982,          4.52246825860982, 3.93437871118446, 6.14896197954192, 6.14896197954192,          3.50829974380329, 3.50829974380329, 3.50829974380329, 3.59391587315661,          3.59391587315661, 3.33646831235876, 3.67933199647321, 3.50829974380329,          3.50829974380329, 4.52246825860982, 4.52246825860982, 3.76454793766266,          3.76454793766266, 4.27163348349591, 4.77149984958079, 4.77149984958079,          4.27163348349591, 4.27163348349591, 3.50829974380329, 3.50829974380329,          3.93437871118446, 3.93437871118446, 3.25025350300496, 3.25025350300496,          3.16383967274129, 4.35544561188301, 4.35544561188301, 3.76454793766266,          4.35544561188301, 4.35544561188301, 5.26417374160422, 5.26417374160422,          3.16383967274129, 3.16383967274129, 3.16383967274129, 4.43905717652239,          4.43905717652239, 4.35544561188301, 4.35544561188301, 3.50829974380329,          3.50829974380329, 3.50829974380329, 4.60567896791123, 4.60567896791123,          4.27163348349591, 4.27163348349591, 4.68868944268447, 4.68868944268447,          4.68868944268447, 4.68868944268447, 4.10340735387365, 2.81620086292774,          2.81620086292774, 5.34559069482765, 5.34559069482765, 3.42248381273475,          3.42248381273475, 3.76454793766266, 3.76454793766266, 4.85411038352648,          4.85411038352648, 4.77149984958079, 4.77149984958079, 3.59391587315661,          3.59391587315661, 5.10074511800645, 3.76454793766266, 3.42248381273475,          3.84956354899072, 3.84956354899072, 2.55340019612735, 2.55340019612735,          3.59391587315661, 3.59391587315661, 3.93437871118446, 3.93437871118446,          4.35544561188301, 4.35544561188301, 3.42248381273475, 3.42248381273475,          4.01899333351728, 4.01899333351728, 4.18762073879345, 4.18762073879345,          5.42680991723477, 5.42680991723477, 3.76454793766266, 3.76454793766266,          2.81620086292774, 2.81620086292774, 6.85553891121225, 6.85553891121225,          5.0187327527967, 5.0187327527967, 3.59391587315661, 3.59391587315661         ), upr = c(14.4208737134878, 14.4208737134878, 13.2880827203392,          13.2880827203392, 11.6111556534885, 13.9812403915557, 13.9812403915557,          12.5235467500559, 11.7752467916346, 15.1343548219701, 15.1343548219701,          11.3665083107223, 11.3665083107223, 12.3558556149815, 12.3558556149815,          12.3558556149815, 12.3558556149815, 14.4208737134878, 12.3558556149815,          12.3558556149815, 12.4396008815305, 12.4396008815305, 12.9462852897287,          12.9462852897287, 11.9401351379537, 11.9401351379537, 12.6076932540179,          12.6920403981972, 12.6920403981972, 12.2723108882934, 12.2723108882934,          12.8613364829343, 12.8613364829343, 12.4396008815305, 12.4396008815305,          11.9401351379537, 12.7765881586919, 11.9401351379537, 11.9401351379537,          13.4601801288326, 13.4601801288326, 12.3558556149815, 12.3558556149815,          12.3558556149815, 12.3558556149815, 12.1058226630474, 12.1058226630474,          12.4396008815305, 12.8613364829343, 12.8613364829343, 12.1889666107397,          12.1889666107397, 11.5294083001765, 11.5294083001765, 12.0228788974822,          12.0228788974822, 12.1058226630474, 12.1058226630474, 12.6076932540179,          12.6076932540179, 13.9812403915557, 13.9812403915557, 12.9462852897287,          12.9462852897287, 11.2043986770552, 11.2043986770552, 12.1058226630474,          12.1058226630474, 12.0228788974822, 12.0228788974822, 11.5294083001765,          11.5294083001765, 12.6076932540179, 12.6076932540179, 12.7765881586919,          12.7765881586919, 12.9462852897287, 12.9462852897287, 15.2243975674713,          15.2243975674713, 12.0228788974822, 12.0228788974822, 11.8575911801404,          11.8575911801404, 12.1889666107397, 13.0314344693091, 12.1889666107397,          12.1889666107397, 12.4396008815305, 12.4396008815305, 12.3558556149815,          12.3558556149815, 13.2023333654031, 13.2023333654031, 14.7760802807705,          14.7760802807705, 12.5235467500559, 12.5235467500559, 12.0228788974822,          12.0228788974822, 12.9462852897287, 12.9462852897287, 12.3558556149815,          14.6869893464325, 14.6869893464325, 11.9401351379537, 11.9401351379537,          11.9401351379537, 12.0228788974822, 12.0228788974822, 11.7752467916346,          12.1058226630474, 11.9401351379537, 11.9401351379537, 12.9462852897287,          12.9462852897287, 12.1889666107397, 12.1889666107397, 12.6920403981972,          13.2023333654031, 13.2023333654031, 12.6920403981972, 12.6920403981972,          11.9401351379537, 11.9401351379537, 12.3558556149815, 12.3558556149815,          11.6931017121067, 11.6931017121067, 11.6111556534885, 12.7765881586919,          12.7765881586919, 12.1889666107397, 12.7765881586919, 12.7765881586919,          13.7198188066704, 13.7198188066704, 11.6111556534885, 11.6111556534885,          11.6111556534885, 12.8613364829343, 12.8613364829343, 12.7765881586919,          12.7765881586919, 11.9401351379537, 11.9401351379537, 11.9401351379537,          13.0314344693091, 13.0314344693091, 12.6920403981972, 12.6920403981972,          13.1167838834176, 13.1167838834176, 13.1167838834176, 13.1167838834176,          12.5235467500559, 11.2853549077749, 11.2853549077749, 13.8067617423288,          13.8067617423288, 11.8575911801404, 11.8575911801404, 12.1889666107397,          12.1889666107397, 13.2880827203392, 13.2880827203392, 13.2023333654031,          13.2023333654031, 12.0228788974822, 12.0228788974822, 13.5465276525046,          12.1889666107397, 11.8575911801404, 12.2723108882934, 12.2723108882934,          11.0430759079299, 11.0430759079299, 12.0228788974822, 12.0228788974822,          12.3558556149815, 12.3558556149815, 12.7765881586919, 12.7765881586919,          11.8575911801404, 11.8575911801404, 12.4396008815305, 12.4396008815305,          12.6076932540179, 12.6076932540179, 13.8939024088035, 13.8939024088035,          12.1889666107397, 12.1889666107397, 11.2853549077749, 11.2853549077749,          15.4956514146983, 15.4956514146983, 13.4601801288326, 13.4601801288326,          12.0228788974822, 12.0228788974822)), .Names = c("Chi", "Sex",      "Age", "SBR", "Diagnosis", "fit", "lwr", "upr"), row.names = c(NA,      201L), class = "data.frame")      I plot SBR v Age for each Sex using `ggplot2`               p <- ggplot(sbr_with_pred, aes(x=Age, y=SBR)) + geom_point(aes(col=Sex),                                                                 shape=19, alpha=0.4) +                  geom_smooth(aes(col=Sex),method = 'lm', se=FALSE,linetype=2) +                  geom_ribbon(aes(y = fit, ymin = lwr, ymax = upr, fill = 'prediction'),                              linetype =2,alpha = 0.1) +                  scale_fill_manual('Interval', values = c('blue')) + theme_bw() +                  theme(legend.position = "right") +                  scale_y_continuous(limits = c(-3,15.5),breaks = c(0,5,10,15)) +                  scale_color_manual("Sex", values = c('red','blue'))      which gives the following ![enter image description here](http://i.stack.imgur.com/j3DxP.png) I can get the equation of each regression fit easy enough               lm(formula = SBR ~ Age, data = subset(sbr_with_pred, Sex == "F"))     lm(formula = SBR ~ Age, data = subset(sbr_with_pred, Sex == "M"))      However how do I test whether or not they are significantly different (which they are not). I think analysis of covariance is the appropriate test but I do not know how to implement this in R


title: Confusion regarding weighted average

text: This is a pretty simple question and I'm wondering how to go about finding a solution. I have the following data.               ttt = data.frame(name=c("A","B","C","D","E","F"),                      count=c(150,250,350,550,150,50),                      returns=c(10,50,60,100,80,10),                      calls_bad=c(5,30,20,15,15,20),                      weight=c(0.20,0.30,0.40,0.40,0.20,0.10))     ttt      For each name, I'm trying to find the return rate, which is just returns/count. However, I want to know how to weight that by how high or low count is. Name F, which only has a count of 50, should be weighted less (weight is 0.10) than name D, who has a count of 550 (weight is 0.40).               ttt$returns_rate = round(ttt$returns/ttt$count, 2)     ttt           name count returns weight returns_rate         1    A   150      10    0.2         0.07         2    B   250      50    0.3         0.20         3    C   350      60    0.4         0.17         4    D   550     100    0.4         0.18         5    E   150      80    0.2         0.53         6    F    50      10    0.1         0.20     d1 = ddply(ttt, .(name), function(x) data.frame(score=weighted.mean(x$returns_rate, x$weight)))     d1       name score     1    A  0.07     2    B  0.20     3    C  0.17     4    D  0.18     5    E  0.53     6    F  0.20      Here are the weights:               Count <= 100 = 10%     Count >= 101 | Count <= 200 = 20%     Count >= 201 | Count <= 300 = 30%     Count >= 301 = 40%      Given what I'm trying to achieve, I thought I'd want a weight average but end up with the following results. However, this can't be right as the returns_rate is equal to the score. How can I go about getting the rate while accounting for the number of count? Thanks! I must've fallen asleep in basic math on this one. EDIT: In lieu of the comments = Eventually, the plan is to combine the call_bad rate (calls_bad/count) and the returns rate (returns/count) into one 'score'/'metric' for each name. However, because count varies significantly by name, I'm working under the assumption that I need to 'weight' the data in order to account for the impact of small 'counts'. Basically, I want to find (returns/count) and (calls_bad/count), and combine these values into one value, while accounting for the fact that someone with a count of 50 could influence the data in a bad way, thus why I'm thinking of using weights. EDIT 2:               ttt = data.frame(name=c("A","B","C","D","E","F"),                      count=c(150,250,350,550,150,50),                      returns=c(10,50,60,100,80,10),                      calls_bad=c(5,30,20,15,15,20),                      weight=c(0.20,0.30,0.40,0.40,0.20,0.10))          ttt$returns_rate = round(ttt$returns/ttt$count, 2)         ttt$calls_bad_rate = round(ttt$calls_bad/ttt$count, 2)          ttt      There only two numbers to "combine" or "average", 'returns rate' and 'calls bad rate'.               ttt$combined = round(ttt$returns_rate + ttt$calls_bad_rate / 2, 2)      But given how count varies "significantly" by name, I thoughts weight were appropriate based on the number of count.


title: How do you get p-values for predictors in a Cox regression analysis in R?

text: I am running a Cox regression analysis in R where both my predictors are binary (`Animal_Number`: _single_ vs. _paired_ ; and `Treatment`: _control_ vs. _exposed_ ). I therefore have four groups ( _single control_ , _paired control_ , _single exposed_ and _paired exposed_ ). When I run my Cox regression, I get the following output (excerpt):                                                       coef exp(coef) se(coef)      z Pr(>|z|)        Animal_NumberPaired                   0.2617    1.2991   0.3303  0.792  0.42829        TreatmentExposed                     -1.4031    0.2458   0.4715 -2.976  0.00292 **     Animal_NumberPaired:TreatmentExposed  1.0246    2.7861   0.5817  1.762  0.07814 .      As I understand it, the three lines presented here are simply the levels of the factors being compared to the baseline experimental group, which was set as _single control_ , i.e., the first line is comparing single control animals with paired control animals. If this is the case, how then does one find information (`coefs`, `z`, `p-value`, etc.) for purely the effect of `Animal_Number`, `Treatment` and their interaction? I can seem to do this in SPSS, and get the following:                                          B    SE      Wald    df  Sig.   Exp(B)     Animal_Number             .244  .331    .545    1   .461   1.276     Treatment               -1.283  .472    7.399   1   .007    .277     Animal_Number*Treatment   .959  .582    2.717   1   .099   2.609      But I want to be able to do it in R. I considered using likelihood ratio tests, to compare full and reduced models, but then information about the coefficients and z is lost.


title: Significant difference between AIC in generalized mixed models

text: I have fit several generalized mixed models for a multinomial distribution data in SPSS. I want to select the best one for explaining my data. But comparing the AIC and selecting the lower AIC model gives me a model without sense and its accuracy is lower than others with higher AIC. I want to test if there is any statistical difference between the AIC of the different models, because maybe the values are not statistically different and I could select another model with higher accuracy. Do you know how to do that? E.g.: `AICa= 1193 (accuracy 81,3%), AICb= 1273 (88,5%)` I know that with R you can test the differences between models, but I can't develop a generalized mixed model with multinomial function in R. Is that possible in SPSS?


title: Floating point error while computing correlation for transformed data using proc corr

text: I get a `floating point overflow` error while using `proc corr` to compute correlations among transformations (exponential and log) of variables. How to avoid this?


title: Statistical properties of a 3D field from spatial averages at different scales

text: I'm trying to estimate the statistical properties of a 3D vector field (the magnetic field vector in the solar atmosphere) over a 2D field of view from observations. I will neglect for the moment that the problem is complex because I cannot measure directly the vector field but only some quantities (the Stokes profiles in some spectral lines) that are non-linearly proportional to the components along the line of sight and on the perpendicular plane. So, I will assume that I can measure $B_z$ (the component along the vertical) and $B_x$ and $B_y$ (in the perpendicular plane). The main difficulty is that the spatial resolution of our best instruments integrates all the information in a square of side $L$, so I get the average of the components on the resolution element $\langle B_z \rangle_L$, $\langle B_x \rangle_L$ and $\langle B_y \rangle_L$. My question is if there is any way in which I can relate the properties of such averages at different sizes $nL$ (multiples of the smallest resolution element) with the statistical properties of the vector field at full resolution?


title: How to 'intelligently' bin a collection of sorted data?

text: I am trying to intelligently bin a sorted collection. I have a collection of $n$ pieces of data. But I know that this data fits into $m$ unequally sized bins. I don't know how to intelligently choose the endpoints to properly fit the data. for example: Say I have 12 items in my collection, and I know the data will fit into 3 bins:               Index:  1 2 3 4 5 6 7 8 9 10 11 12     Value:  1 1 1 3 3 3 3 3 3 5  5  6      How do I intelligently choose my breakpoints for the bins of $i = \\{1-3\\}, \\{4-9\\}, \\{10-12\\}$? The current implementation I have breaks the data into evenly sized bins and then takes the average of the endpoints to find the indexes for the end of the bins. So it works like this:               Index:  1 2 3 4 5 6 7 8 9 10 11 12     Value:  1 1 1 3 3 3 3 3 3 5  5  6          first break evenly: i = 1-4, 5-8, 9-12     mean endpoints:  between 4 and 5: (3+3)/2 = 3                      between 8 and 9: (3+3)/2 = 3      So now anything below 3 fits in bin 1, anything above 3 but below 3 fits in bin 2, and anything above 3 fits in bin 3. You can see what my problem is. If the data has unequal bins my method fails. A friend mentioned the k-nearest neighbor algorithm but I'm not sure.


title: How to check psychometric type of scales?

text: It is well known, that there are diffrent types of scales (see Wikipedia Level of measurement). In psycholinguistics there is Semantic Differential technique. It uses a scale (usually 5 or 7 grades) which type as **ordinal** or **interval** is being debated. As Wikipedia says, the majority think it is ordinal. But at the same time processing researchers use arithmetic mean and Factor Analysis: and it uses `+` and `-` math operations. Is it critical error? May be, nevertheless, in fact the intervals between closest grades are equal - how to check it? What is the cause of scale's type?


title: According to the Neyman-Pearson framework is it incorrect to say something 'trends towards significance'?

text: I've had it put to me by a couple of statistics teachers that statements of this kind are just loose talk. This seems to relate to the discussion of the Neyman-Pearson framework I've read about at When to use Fisher and Neyman- Pearson framework?. Did Neyman and/or Pearson ever criticize the use of statements like "this result trends towards significance"? If not, are such criticisms of statements like this nonetheless inherent to the Neyman-Pearson approach? Did Fisher ever endorse the use of such statements?


title: Uniformly Most Powerful (UMP) test for Uniform(t,t+1)

text: Suppose $X_1,...,X_n$ follow a $U(t,t+1)$ distribution. Show that the UMP level alpha test for testing $H_0:t<=t^*$ vs $H_1:t>t^*$ has form: Reject $H_0$ if and only if $\min(X_1,...,X_n)>t^* + c(\alpha)$ or $\max(X_1,...,X_n)>t^*+1$ for some appropriate $c(\alpha)$.


title: How to validate a Bayesian result?

text: I'm trying to link statistical analysis methods to validable business observations. Ideally an interpretation that is measureable and can objectively rank success of predictions. For the frequentist p-value of 0.95 I figured the statement is: There are instances - some of which have a real effect and some are purely random (according to some agreed notion of randomness). Amongst the purely random instances, this frequentist will misclassify 5% as having an effect. (Nothing else can be said?) I think I can assume that p-value statements transfer to a set of instances, even if it's completely different experiments(?). No need for repeatable experiments. Is this correct so far? What's an equivalent statement for Bayesian results? It should be something measureable.


title: Comparing and contrasting classification trees over the years

text: At work I have decided to use a classification tree to do some analysis on the people at the bottom decile for hourly pay. Mainly for inferential reasons. My datasets are broken up in years, I have 2013, 2012, ... , 2006. I just wanted to know if I can build trees for each of those years, and then compare and contrast them? So in one year the total hours worked variable may be in there as one of the nodes, and in another year it may not be, can I conclude that in the year that it was, it was a more important factor in the person being in the bottom decile? I'm using the same methodology for them all, i.e the datasets for each year are set up the same way with the same variables, no missing values, cross- validating to get the best cp value, and then building the tree. Thanks in advance. Rafi


title: Binomial regression and heteroscedasticity-consistent (HC) standard errors

text: I have a small (n = 60) dataset where I would like to estimate the magnitude of the relationship between a set of predictor variables (among others) and an outcome rate variable. Because the number of trials for the outcome variable differs dramatically between cases, it seems natural to weight each case by the number of trials. Perhaps more problematically, for theoretical reasons we expect that the number of cases is correlated with some of our predictor variables. It seems like a natural fit to use a binomial regression here. However, I am worried that heteroskedasicity may cause the binomial regression to underestimate the standard errors. **Should I report heteroscedasticity- consistent (HC) standard errors here? If so, can I use the`sandwich` package in R to estimate them, or will this mess with the weights specified in the `glm` call to account for the differential number of trials?**


Does it make sense to compute confidence intervals and to test hypotheses when data from whole population is available? In my opinion, the answer is no, since we can accurately compute true values of paramaters. But then, what is the maximal proportion of data from original population that allows us to use aforementioned techniques?


title: Statistical test for difference between two odds ratios?

text: I am attempting to compare two diagnostic odds ratios (DORs). I would like to know of a statistical test which will allow me to do this. Please help! Thank you!


title: Is R-square enough to show that my data fitted well to my regression tree model?

text: In the toy example given below, is it enough to show that the model fitted weel? Is there any other metric available to show the fitness? Thanks               library(rpart)     N <- 1000     X1 <- runif(N)     X2 <- 2*runif(N)     X3 <- ordered(sample(letters[1:4],N,replace=TRUE),levels=letters[4:1])     X4 <- factor(sample(letters[1:6],N,replace=TRUE))     X5 <- factor(sample(letters[1:3],N,replace=TRUE))     X6 <- 3*runif(N)      mu <- c(-1,0,1,2)[as.numeric(X3)]     SNR <- 10 # signal-to-noise ratio     Y <- X1**1.5 + 2 * (X2**.5) + mu     sigma <- sqrt(var(Y)/SNR)     Y <- Y + rnorm(N,0,sigma)     # introduce some missing values     X1[sample(1:N,size=500)] <- NA     X4[sample(1:N,size=300)] <- NA     data <- data.frame(Y=Y,X1=X1,X2=X2,X3=X3,X4=X4,X5=X5,X6=X6)          #Divide into training and test data      testidx = sort(sample(1:N, 200, replace=FALSE))     dataTesting <- data[testidx,]     dataTraining <- data[-testidx,]     #Fit the data using Rpart     fit <- rpart(Y~X1+X2+X3+X4+X5+X6, method="anova", data=dataTraining)          #prediction     predicted <- predict(fit, newdata=dataTesting[,-1], type="vector")          #Find R-square     actual <- dataTesting$Y     rsq <- 1-sum((actual-predicted)^2)/sum((actual-mean(actual))^2)     rsq


title: How to perform RMSE analysis in SPSS?

text: My thesis coach wants me to perform a predictive analysis based on OLS. What I understand is the following:   * divide the dataset into a training set and a holdout set, for instance 50-50   * perform OLS on the training set   * construct linear equation based on regression output   * create a new variable (DV2) in the holdout set, and use the linear equation to calculate its values   * now you have F (forecasted) and A (actual) DV values in the holdout set   * calculate the performance of the predictive linear equation with RMSE: ![enter image description here](http://i.stack.imgur.com/YS8le.png)   * a lower RMSE is better Now, some questions about this procedure:   * Am I doing this right?   * I have no clue how to have SPSS perform the RSME operation, so can't I just do it in Excel? If I paste the holdout set into Excel, performing this calculation seems easy enough. Is there something I'm missing?   * If you know how to perform this calculation in SPSS, please let me know because I expect that SPSS might be able to output some extra insightful statistics and / or graphs


title: What terms should I include in a linear regression model?

text: Consider the following OLS output. My question is: If you have to define a linear equation from the regression, do you then only include significant variables? If so, the equation according to the SPSS output (DV = `stockVol0LN`) would be:               stockVol0LN = 0.632 * StockVol1LN + 0.278 * StockVol2LN + 0.152 * hbVol0LN -                    0.050 * hbVol1LN - 0.069 * hbVol2LN - 0.045 * hbAgreeQuality0LN +                   0.064 * wiki0LN - 0.067 * wiki1LN + 0.050 * svi0 + 1.396      ![enter image description here](http://i.stack.imgur.com/FTGPk.png) **My question is:** Is this the best equation, or is there another measure besides significance level to consider when deciding what terms you should include in the model?


title: Doubts on principal component analysis

text: I need to reduce the dimensionality of motion capture data $M=\\{x,y,z,\theta\\}$ where $x,y,z$ are the pixel coordinates and $\theta$ is the Euler angle/joint angle. This data is obtained from skeletal markings in video frames. So, for 10 joints I have 10 such sets of time series data giving a total of 4*10=40 time series data. I need to apply PCA to reduce the dimension to maybe say 1D time series data. The doubt are (A) Do I need to convert/transform the data from video coordinate to world coordinate system? (B) after doing PCA on $M$, the new representation on which I will do further processing will be the eigenvalues or the original raw data after the main component has been identified which contributes maximum in the data set? How do I apply PCA in this case? A simple illustration of how PCA will work will be really useful.


title: In regression analysis, if a potential third variable is a function of two others, is it (necessarily) redundant?

text: Say I am using regression analysis to try and estimate the assembly time for widgets (which can be produced in various sizes) given some possible independent variables:   * Width   * Height   * Quantity of pieces (different sizes used as components)   * Did the American League win the World Series   * etc. Most finished widgets are rectangular, however I'm wondering if the "squareness" of a widget explains some of the assembly time. That is, let's say Widget A is 2 feet wide and 18 feet long. Widget B is 6 feet square. Each have the same area (36 sq ft), but my hunch is that because of the assembly process, Widget B will take longer due to its square shape. Now to my question... Since "squareness" is really just a calculated ratio of two other included variables (width and height), is it already "explained" in these two variables? Is it redundant to include "squareness" or similar ratio as an independent variable?


title: overall effects of categorical variables

text: I'm doing a Poisson regression in Stata, so the dependent variable is a count variable and I have some categorical predictors. If A is a categorical variable with, for example, 4 levels, in the parameter estimates table I get results for the 3 levels of the variable compared to the level I have set as the reference category. Similarly for an interaction term. However, I would also like to have an estimate of the overall effect of variable A and also an estimate of the overall effect of the interaction. How do I do this in Stata? Ie. I need the test of the effect of A in the model on the dependent variable D, which tests the joint hypothesis that every one of the four levels of A has the same effect on D, and therefore A does not explain any of the variation in D, and so it is not a significant predictor of D in the model. I do not need a test which tests the hypothesis that the effects on D of e.g. A2 , A3, or A4, respectively, are all equal to each other in this model, but not necessarily equal to the effect of A1. Thank you UPDATE: Well, prompted by the comments from all of you, and after hours of trying out different things, I found that what I was looking for was "constrast A", described in http://www.stata.com/manuals13/u25.pdf However, I am now wondering why "contrast A" does not give me the same results as "test A4=A3=A2=A1". I'm told I shouldn't make stata-specific questions, however I'm trying to see how stata corresponds to other packages and this is what prompted my question. I also think that if someone can answer this question for me i.e. why constrast is not the same as a joint test on whether all levels have the same effect, then this would be of general interest as it obviously translates to a general statistical question.


title: How to test for randomness in bins with small N?

text: I observe a series of crime incidents linked by _modus operandi_ or some other peculiar characteristic of the crime (e.g. cutting catalytic converters from underneath vehicles). I would like to know if the observed days of the week that the crimes occur on (ignoring the uncertainty that sometimes occur in crime incidents - e.g. it happens overnight) are random. I typically have very few linked crime events, say between 5-15. So question 1, there is a lot of knowledge about using Pearson's $\chi^2$ on small n contingency tables, can I use all of that same advice for $2\times 2$ contingency tables that come with it for just 7 day of week bins? (In particular can I use the $N - 1$ correction given expected cell frequencies are at least 1 and still expect similar coverage rates - which would mean I need at least 7 observed events?) Or alternatively question 2, are there any other obvious approaches I can take to test the hypothesis of the events being random with respect to the day of the week? (Permutation approaches perhaps given the limited nature of the potential permutations?)


title: Does stand-alone dummy variables in linear regression models make sense?

text: Dummy (or binary) variables ($X_2$) can be used in linear regression models to help explaining a possible group effect that a continuous predictor variable ($X_1$) might present in explaining the response variable ($Y$). Now, I am wondering if it makes sense to have only dummy variables as predictors in a linear regression model with a continuous response variable. Does it? Example: % of population with instruction = influenced by politics A + influenced by politics B. Both politics A and B can assume values 1 or 0.


title: Is non-stationary AR(p) process constant in mean?

text: A non-stationary $AR(1)$ process, which is a random walk, is constant in mean, but not constant in variance. How about the other $AR(p)$ processes with the order $p>1$? Are they constant in mean?


title: Hyperprior distributions for the parameters (scale matrix and degrees of freedom) of a wishart prior to an inverse covariance matrix

text: I'm estimating several inverse covariance matrices of a set of measurements across different subpopulations using an wishart prior in jags/rjags/R. Instead of specifying a scale matrix and degrees of freedom on the inverse covariance matrix prior (the wishart distribution), I would like to use a hyperprior on the scale matrix and degrees of freedom, so they can be estimated from the variation across subpopulations. I haven't found much literature on hyperpriors for the scale matrix and degrees of freedom. Most of the literature seems to stop the hierarchy at the choice of the prior to the covariance/inverse covariance and/or are focused on estimating a single covariance matrix rather than several covariance matrices across different populations. Any suggestions as to how to go about this - what are the recommended hyperprior distributions to use for the scale matrix and degrees of freedom of the wishart distribution? Is there some literature on this that I'm missing?


title: Statistical significance of breakpoints from `segmented` in R

text: I'm looking for breakpoints in species abundance as a function of spatial distance, using the `segmented` package for R. 'segmented' appears to return a breakpoint no matter what; I don't understand whether it returns an estimate of the signifcance of the breakpoint (or whether the segmented linear model is better than an unsegmented model). For instance (R code):               require(segmented)     set.seed(1)     x <- 1:100     y <- rnorm(100) + x # No real breakpoint     y2 <- c(50-x[1:50], x[51:100-50]) + rnorm(100) # Clear breakpoint at 50     plot(x,y)     points(x, y2)      ![enter image description here](http://i.stack.imgur.com/DRysM.png)               # Segmented model for the unsegmented data     testM <- lm(y ~ x)     testMs <- segmented(testM, seg.Z = ~x, psi=90)     summary(testMs)      Despite the fact that there is clearly no breakpoint in the data, `summary` reports that `t value for the gap-variable(s) V: 0`, and the standard error is fairly small (3.28).               testM2 <- lm(y2 ~ x)     testM2s <- segmented(testM2, seg.Z=~x, psi=50)     summary(testM2s)      Here the guess is correct (49.78), the standard error is even smaller (0.2), and the t value is the same. How do I interpret this result? **Note:** I have no particular attachment to the `segmented` package - I just want to test the hypothesis that a segmented regression models my data better than a single-domain regression. But the other breakpoint-analysis packages I've looked at seem to require that points on the domain be evenly-spaced (e.g. timeseries) and there is a single dependent value per independent value. These assumptions are not met for my spatial data.


title: MLE of a Cauchy distribution

text: > After centring, the two measurements x and −x can be assumed to be > independent observations from a Cauchy distribution with probability density > function: > > $f(x :\theta) = $ $1\over\pi (1+(x-\theta)^2) $ $, -∞ < x < ∞$ > > Show that if $x^2≤ 1$ the MLE of $\theta$ is 0, but if $x^2>1$ there are two > MLE's of $\theta$, equal to ±$\sqrt {x^2-1}$ I think to find the MLE I have to differentiate the log likelihood: $dl\over d\theta$ $=\sum $$2(x_i-\theta)\over 1+(x_i-\theta)^2 $ $=$ $2(-x-\theta)\over 1+(-x-\theta)^2 $ + $2(x-\theta)\over 1+(x-\theta)^2 $ $=0$ So, $2(x-\theta)\over 1+(x-\theta)^2 $ $=$ $2(x+\theta)\over 1+(x-\theta)^2 $ which I then simplified down to $5x^2 = 3\theta^2+2\theta x+3$ Now I've hit a wall. I've probably gone wrong at some point, but either way I'm not sure how to answer the question. Can anyone help?


title: Are the predictors in an ANOVA assumed to be independent of one another?

text: I ran an ANOVA with 4 predictors (2 categorical, 2 continuous) and 1 covariate (continuous). Including the covariate did not significantly alter the main effects of the ANOVA when it was run without the covariate. I then found that one of the categorical variables had a significant effect on my covariate (i.e. the covariate was significantly higher at one level of the categorical variable compared to the other). 1) Does this mean that these predictors were not independent? 2) For that sort of ANOVA, was I meant to assume that all the predictors would be independent? 3) What are the consequences for the interpretation of the main effects if these variables were not independent?


title: Using the predicted values of one regression as the dependent variable in another; problems? techniques?

text: Forgive me if this question is redundant; predicted values and regression returns many results. My question is this: I would like to run one regression and generate predicted values based upon a subset of the variables in that regression. With these predicted values I then want to regress them on other variables. Multi-level modelling is similar to the task, but I specifically need the predicted values as opposed to the coefficients. Based upon this strategy, would I need to do anything with the error structure in the second regression? FGLS with a het correction based upon the subset of variables seems reasonable. But i'm out of my depths in terms of theory. Tia.


title: Computing (lagged) correlations (or similar) between multiple time-series from a VECM or its levels-VARM

text: Ok, after trying to find a way to get all lagged and non-lagged - correlations from multiple time-series (behavioral series with length around 180) I've confined myself to the following: I have multiple (5 so far, more coming) sets of 12 time-series with equidistant measurements, same sampling-times $t_1, t_2, ... ,t_{180}$. Not all of them are stationary, but I'm assuming for the moment, most if not all of those will be participating in some cointegration-relation. (Will have to check that, but just assuming for the time being...) For each set (containg 12 parallel time-series) I would like to estimate a VECM and, based on that, would like to get some unique measure of intra- and inter-series correlations for lags=0..p (where p is the order of the VECM - i expect p to turn out 2..3, 4 at the max). Now, all I'm really interested in are **unique** parameters of "connection", it doesn't have to be actual correlation, as long as it is unique and thus comparable between different sets of the 12 series. (I would like to compare a certain subset of the parameters of one of the sets of series with the parameters of the other sets of series). For now, I also don't really care, wheter the parameters are partial (correlations) or not. My idea is to estimate the VECM and (supposing I can get a good fit for the model by some criterion) reformulate the VEVM(p) as a levels-VAR(p+1) and then using the VAR's Matrix-Parameters as measures of (lagged) "connections" between the series. Only problem (I see) is: How can I ensure uniqueness of the estimated VAR-model's parameters - I know the cointegration vectors in the VECM are unique only wrt. scalar multiplication by some constant, so I could just normalize the cointegration vectors to length=1. Would this yield a unique VAR? Also, I've found this post, which gave me the idea to Z-transform the series beforehand, I'm quite new to the cointegration/VECM/VAR-business, but I would assume this (combined with normalized cointegration vectors) should give me unique VAR-parameters. Does anyone have some advice on this, confirmation or criticism?


title: Interpreting interactions in multilevel model

text: I'm running a rather simple multilevel model in which I have two variables interacting at level 1. The interaction term is significant, but I'm not sure how to interpret it properly. I am familiar with doing interaction analyses in linear regression, do I do the same steps? I'm using Stata, in case that's of any help, or there are any cool tools (such as sslope) for making this easier.


title: What are the most popular references for mixed models?

text: What are the most popular references (especially texts) you think in the area of mixed models? Or what are those you think were helpful to you when you learnt this topic?


title: Support as a Proportion in Apriori

text: ![enter image description here](http://i.stack.imgur.com/rrHOF.png) I'm asked to find association rules with support >= 0.3 and confidence >= 0.6. For support as a percentage, we need to divide the support of the candidate itemset by the the total number of records at each k-item stage of the algorithm? ie., for 1 item itemsets: supp({1}) = 3/8 supp({2}) = 5/8 supp({3}) = 5/8 all pass the threshold, so, for 2 item itemsets: supp({1,2}) = 2/8 supp({2,3}) = 2/8 supp({1,3}) = 3/8 Giving 1 -> 3 and 3 -> 1 as possible association rules? (need to check confidence).


title: Appropriateness of applying a fitted model to a different (but similar) set of predictor variables

text: I have fitted a model of current species habitat suitability as a function of mean annual temperature (MAT) and precipitation (MAP) by regressing the distribution of known occurrences of a species with respect to these variables across space using GBM. Now I have a colleague that wants to employ this model to look at seasonal shifts in habitat suitability by applying the model to mean temperature and precipitation data for summer and winter months. In my mind, this is not a valid use of the model as both climate and species location data was pooled across seasons when fitting the model. As I understand it, swapping them does not answer a question about the distribution of a species in the dry season, it instead answers the question of what the distribution of the species would be if MAT was as high as dry season MT. My question is, could this type of substitution be ever appropriate, especially when the model projection variable is just a subset of the original predictor variable? If not, what are the principles/concepts that most easily explain why this should not be done?


title: Test whether two rank orders differ

text: I have 20 human subjects doing a task under 2 different difficulty conditions. (So it's a repeated measures design) Each person has a different performance in each of the task difficulty conditions, so I can rank-order them. I want to test whether there are differences between the ranking of performaces between conditions. This question might be unbelievably stupid, but I'm essentially looking for a test that would give me a significant result when there is no correlation between measures. EDIT: To clarify, I want to see if the rank order 1 2 3 4 5 6 7 8 9 10 is different from the rank order 1 3 4 5 2 8 9 6 7 10 The Wilcoxon signed-rank test -and other rank-based statistics- won't do, because their null hypothesis (unless I'm mistaken) is that the median of the two ranks is equal. This allows to test overall differences in rank order (when the two conditions are ranked together). What I need is a test that will compare whether the order itself changes. Can anybody point me to a test that has the null hipothesis that the rank _orders_ are equal? Thanks


title: Mechanism of multiple imputation?

text: I am trying to understand the mechanism underline multiple imputation ideas. I am confusing on creating multiple sets of estimates and then average them. For example:               set.seed(123)     df<-data.frame(     x1=rnorm(100,10,2),     x2=rbinom(100,20,0.6))          y=rbinom(100,1,0.5)     #set missing values     ms<-sample(1:100,20)     y[ms]<-NA     df$y<-y          #predict NA          fit<-glm(y~x1+x2,data=df,family="binomial")     pred<-predict(fit,newdata=df[,c("x1","x2")],type="response")     df$yfit1<-with(df,ifelse(pred>=0.4,1,0))         df$yimp1<-with(df,ifelse(is.na(y),yfit1,y))                   x1 x2  y yfit1 yimp1     1  8.879049 14  1     1     1     2  9.539645  8  0     1     0     3 13.117417 11 NA     1     1     4 10.141017 12  1     1     1     5 10.258575 13  1     1     1     6 13.430130  9  0     1     0      If I want to impute `m` times to find the average estimates for `x1` and `x2`, in what idea to create multiple `df$yimp1...df$yimpm`. Why are not these estimates same in different runs? Could someone give an example for `m=3`?


title: Detecting statistically significant clustering of continuous values

text: I'm working with biological sequence data where each position in the sequence has an associated continuous value. I'm ignoring the sequence content so the data is very similar to a time series with measurements at discrete timepoints -- all values are equally spaced. I would like to be able to detect whether high values tend to cluster together (occur in runs) so I applied the Wald- Wolfowitz runs test for non-random placement of values >1. There are some issues with that approach:   * Wald-Wolfowitz works on binary data so I have to binarise the continuous values I have (everything larger than 1 becomes 1 and the rest is 0). Ideally I would like to be able to detect features such as runs of similar values (let's say 10 values of 0.5 in a row) as well. I would imagine there are some methods that would operate on continuous values (e.g. based on autocorrelation) but couldn't find any.   * While I get a measure of clustering (the test p-value), I don't know which parts are actually clustered or how many clusters there are.    * I would also like to extend this approach to 3D (mapping of sites on the protein structure) and the test doesn't support multiple dimensions, either. I was wondering if there are more sophisticated statistical approaches that I could apply?


title: E-mail answered probability after n days of waiting for a reply - based on a sample of e-mails and replies

text: Here is the task: I have a sample of replies to my e-mail from my mail box. A sample is taken over a period of 90 days, 1000 e-mails and replies if any. (We only consider a pair of {my original e-mail, date}, {reply, date}). Question: how many days should I wait to conclude I will never get a reply from this person (given an acceptable error I am OK with, like in 10% cases there will be a mistake), I can make this conclusion earlier than 90 days waiting period. So, I take my 1000 sample and arrange as follows: number of people that replied within 0 days (on the same day) - x0; replied in 1 day - x1; ... replied in 60 days - x60; ... replied in 90 days - x90; did not reply in 90 days - x91. Now, I can see then, for example, 50% - 500 people replied within 10 days. 75% replied in 15 days; 95% replied in 30 days, etc. Based on my numbers I need to set a reasonable threshold for %. As if in 15 days it was 75% but in 27 days it was 80% I can safely choose 15 days of waiting. Maybe after that many days of waiting I will stop check if that person replied 5 times a day and will do it every other day, for example, or whatever action I want to take based on my observation. How do I formulate this task correctly in math/statistical/probability terms and what's the solution? Do I need to estimate a distribution here? Thanks!!!


title: How to work with conditional expectations and variances?

text: I have seen this statement in a lecture note : $$E[(Y-f(X))^2\,|\,X]=Var[Y|X]+E^2[(Y-f(X)\,|\,X]]$$ If $Z,X$ are random variables, then shouldn't be $Var[Z|X]=E[Z^2|X]-E^2[Z|X]$?In which case shouldn't the term $Var[Y|X]$ be $Var[Y-f(X)\,|\,X]$?


title: Fisher's method when p-value = 0

text: I have a bunch of independent _p-values_ and now I want to combine them using the Fisher's method. However, in R, when the p-value is zero, the `log(p-value)` becomes $-\infty$, so that the test statistic $X^2$ is $\infty$. This will give the Fisher's method _p-value_ of 0. If we check the _p-values'_ histogram, they look pretty uniform (under the null) so that I expect the Fisher's method should give me a not-so-small _p- value_ (should not reject the null). Is there a way to get around this issue? How about the Stoffer's method? Will that method automatically take care of the issue when $p = 0$ (or close to 0)? **Update:** I tried the Stouffer's method, but still the issue of $\infty$ is not solved...


title: cross tabulation of two variables with resulting average and SD

text: I am beginner of R software. I am looking for two variables and result cross tabulated but not in the form of frequencies. Could I get help on how to get result as average and standard deviation instead of frequencies while cross tabulating in R. my data looks like this               Treatment   Stage   Chlorophyll     Salt        Green   0.2     Salt        Green   0.3     Salt        Green   0.4     Salt        Pink    0.5     Salt        Pink    0.3     Salt        Pink    0.2     Salt        Red     0.5     Salt        Red     0.6     Salt        Red     0.7     Nitrogen    Green   0.4     Nitrogen    Green   0.6     Nitrogen    Green   0.9     Nitrogen    Pink    0.2     Nitrogen    Pink    0.3     Nitrogen    Pink    0.5     Nitrogen    Red     0.4     Nitrogen    Red     0.2     Nitrogen    Red     0.3     Control     Green   0.5     Control     Green   0.6     Control     Green   0.8     Control     Pink    0.5     Control     Pink    0.4     Control     Pink    0.6     Control     Red     0.2     Control     Red     0.3     Control     Red     0.1      I appreciate help from this community.


title: Bootstrapping with random effects in SPSS

text: I'd like to use bootstrapping in my two-way ANOVA containing two fixed and one random factor. Why is the bootstrapping method not available (greyed out) for models containing a random factor? Thanks!


title: GARCH(1,1) regression in Eviews

text: I'm having a problem in doing a GARCH(1,1) regression. I'm trying to regress gold prices serie on stock returns series as in the following equation in eviews: $$ r = c(1) + c(2) \cdot s + c(3)\cdot d10\cdot s +c(4)\cdot d05\cdot s + c(5)\cdot d01\cdot s , $$ where $r$ is gold prices, $s$ stands for stock prices and $d$ represents dummy variables to capture the 10% 5% and 1% high volatility. My questions are :   1. Should any of the two series be stationary before I do the regression?   2. Should I write my equation in "the mean equation" or in "the variance equation" or both ? I'd really appreciate if someone could tell me the stages I should go through to solve this problem, I'd be deeply indepted. Note that I'm a beginner in this GARCH modelling, the more I read about it the more I get confused.


title: Basic conditional sampling question

text: Let's say that I have a method of sampling from $P(X, \theta \mid Y)$. It seems intuitive that if I drew samples of $(X, \theta)$ and then simply dropped the $\theta$s, this is equivalent to drawing from $P(X \mid Y)$ directly. What is the theorem that says that I can do this?


title: Inferring prior distribution

text: Suppose that we take a sample ($X_1, X_2, ... X_n$) from a distribution where we assume that $X_i $~$ Bin(n_i, p_i)$ and $n_i$ is known for every $i$. We also assume that $p_i$'s are independent and identically distributed, $p_i$ ~ $D$, where $D$ is some unknown distribution. $n_i$ cannot be assumed to be large. My goal is to get a Bayesian estimate (or a probability distribution) for $p_i$. But this requires coming up with a distribution for $D$. One option is to make an empirical distribution that uses frequentist estimates for each $p_i$ (i.e. $p_i = X_i/n_i$). This is a rather intuitive and potentially reasonable idea. Unfortunately, the presence of small $n_i$'s would make the tails heavier than they should be (lots of extreme values close to 0 or 1). **I'm looking for another option that doesn't have the problems of the aforementioned solution.** One possibility I have in mind is to use the following algorithm:   1. Generate prior distribution as explained earlier.   2. Get MAP or EAP estimate for every $p_i$.   3. Generate new empirical prior from the probabilities obtained in 2.   4. Go back to 2 (continue for a set number of steps, or possibly until convergence?) Is this method similar to any method out there? Is it reasonable?


title: Using LASSO from lars (or glmnet) package in R for variable selection

text: Sorry if this question comes across a little basic. I am looking to use LASSO variable selection for a multiple linear regression model in R. I have 15 predictors, one of which is categorical(will that cause a problem?). After setting my $x$ and $y$ I use the following commands:               model=lars(x,y)     coef(model)      My problem is when I use `coef(model)`. This returns a matrix with 15 rows, with one extra predictor added each time. However there is no suggestion as to which model to choose. Have I missed something? Is there a way I can get the lars package to return just one " _best_ " model? There are other posts suggesting using `glmnet` instead but this seems more complicated. An attempt is as follows, using the same $x$ and $y$. Have I missed something here?:               cv=cv.glmnet(x,y)     model=glmnet(x,y,type.gaussian="covariance",lambda=cv$lambda.min)     predict(model, type="coefficients")      The final command returns a list of my variables, the majority with a coefficient although some are =0. Is this the correct choice of the " _best_ " model selected by LASSO? If I then fit a linear model with all my variables which had coefficients `not=0` I get very similar, but slightly different, coefficient estimates. Is there a reason for this difference? Would it be acceptable to refit the linear model with these variables chosen by LASSO and take that as my final model? Otherwise I cannot see any p-values for significance. Have I missed anything? Does               type.gaussian="covariance"       ensure that that `glmnet` uses multiple linear regression? Does the automatic normalisation of the variables affect the coefficients at all? Is there any way to include interaction terms in a LASSO procedure? I am looking to use this procedure more as a demonstration of how LASSO can be used than for any model that will actually be used for any important inference/prediction if that changes anything. Thank you for taking the time to read this. Any general comments on LASSO/lars/glmnet would also be greatly appreciated.


title: Feature selection with SNP datasets

text: I'm working with a Single-nucleotide polymorphism (SNP) dataset with over 2.2 million features and roughly 2000 samples. I wish to do feature selection on this dataset to reduce it down to approximately ~20,000 features so that I can then use them to train different classifiers for a binary classification problem. The 2 classes are in the ratio of 3:1 What are the best feature selection methods in this scenario?


title: Is mRMR feature selection sensitive to imbalanced dataset?

text: I wish to perform feature selection on a dataset using mRMR but my classes are imbalanced roughly 3:1. Should I down sample the majority class or use the whole dataset? If I use the whole dataset will it select features that are biased towards the majority class?


title: How does the interpretation of main effects in a Two-Way ANOVA change depending on whether the interaction effect is significant?

text: As far as I can tell this question has never been asked before. There are several questions that touch on related issues, but as far as I can see none of them have provided a definitive answer to this question. Furthermore, in at least one case the most upvoted and second most upvoted answer have implicitly disagreed with each other. Textbooks often caution against interpreting main effects in Two-Way ANOVA when the interaction effect is significant. Here's a (mild) example of this, from pp562-563 of Hatcher's "Step-by-step Basic Statistics Using SAS: Student Guide": > When you perform a two-way ANOVA, it is possible that you will find that (a) > the interaction term is statistically significant, and (b) one or both of > the main effects are also statistically significant. When you prepare a > report summarizing the results, you will certainly discuss nature of your > significant interaction. But is it also acceptable to discuss and interpret > the main effects that were significant? > > There is some disagreement between statisticians in answering this question. > Some statisticians argue that, if the interaction is significant, you should > not interpret the main effects at all, even if they are significant. Others > take a less extreme approach. They say that it is acceptable to interpret > significant main effects, as long as long as the primary interpretation of > the results focuses on the interaction (if it is significant). Do I assume correctly that the hard-line stance ("you should not interpret the main effects at all") is incorrect, but that the interpretation must nonetheless change if the interaction effect is significant? If so, how does the interpretation change?


title: How to perform this ANOVA?

text: I am interested to know if we have 3 dependent variables, and 2 independent variables, online and face to face, with 2 levels each, as in the example below. The sample size is 30 for each mode of study - 30 participants for online mode and 30 for face to face.   1. What type of ANOVA do we need to conduct?   2. Do we have to analyze each dependent variable at time? For instance, we first conduct an ANOVA for "graduation Time"; then another ANOVA for "Recall" and finally another ANOVA separately for "Job Skills"?   3. Will it be more complicated to conduct a MANOVA, instead of conducting separate ANOVAs? Example: Mode of study (online (6 months, 12 months), face to face (6 months, 12 months)) has an influence on individual’s graduation time, recall of material, and job skills


title: Analyzing relationship between two Likert-items (is it possible?)

text: I am a student and am new to statistical market research. Is it possible to analyze the dependence of answers to one Likert-type item on answers to another Likert item (predictor)? If so, what test can I use? I am doing an analysis of how much attachment/affinity (dependent Likert item) of consumers to a certain brand can be affected by their perception of how much a new brand extension of that certain brand fits the current brand image they have in mind (independent/predictor Likert item). I am thinking a linear regression can be applicable here. But please suggest other ways if it is not. I am very open to any suggestion and advice. Thank you very much. :)


title: Simple measure for comparison of rainfall regularity

text: In meteorology we have the concept of monthly rainfall, which is just the sum of daily rainfall over that month. Now, given this extreme example: ### Situation 1:               First day of month: 210mm     All remaining days: 0mm          Monthly rainfall = 210 + 0 + ... + 0 = 210mm      ### Situation 2:               All days of month: 7mm     Monthly rainfall = 7 * 30 = 210mm      Which statistic can I use to show that Situation #2 has a more 'regular' rainfall distribution? I'm not sure if standard deviation or entropy will give me a good result.


title: How to find trends/correlations in a simple 2D dataset?

text: I have a dataset of approximately 100 objects. For each object, I have made measurements and have a handful of values (typically less than 10). I also know several physical properties of my objects. My dataset is arranged thusly: _(arbitrary numbers)_                         Measurement1  Measurement2  Measurement3 ... Property1  Property2 ...     Object1       42.08         85.23         29.44          8000K      0.001     Object2       27.89         53.26         20.85          6200K      0.008       ...       etc      It is known that some trends are present. For instance, linear regression of _Measure1_ vs. _Measure2_ (when both are normalized to _Measure3_ ) show that these quantities are well correlated. This is reproduced in my data when measured by hand. The idea is to find trends or relationships that are not already known (to be frank, we know little at present). How might I accomplish this? An example would be very appreciated. Python is preferred.


title: Advantage of kernel density estimation over parametric estimation

text: Is there any particular reason you will choose the kernel density estimation over the parametric estimation? I was learning to fit distribution to my data. This question came to me. My data size is relavtively large with 7500 data points. Auto claims. My goal is to fit it to a distribution(nonparametric or parametric). And then use it to simulate auto claim data, and calculate VaR or TVaR. I used log to transform the data to make it relatively normal. I fitted many distributions including normal, lognormal, gamma, t,etc...I used AIC and loglikehood to identify the best fitting. But none of all this fitting passed KS test(p value extremmly small, with e-10). That's why I asked in what situation I should switch to KDE.


title: Support Vector machine : a simple question

text: I think that a formulation of SVM for points x with label y is : $$ \begin{align} \arg\min_{\substack{u,w,b}} \frac{1}{2} \cdot |w|^2 + C \cdot \sum_{i} u_i \\\ s.t.\ \ y_i\cdot (w \cdot x_i + b) &\geq 1-u_i \\\\\\\ u_i \geq 0 \\\ \end{align} $$ In that formulation, if we take C = 0, what prevents $u_i$ to go to infinity, so that the two constraints are always satisfied ?


title: Hyper-parameter optimization via random search

text: I’m working on a classification system which consists of an auto-encoder for feature learning and logistic regression for classification. The system has five hyper-parameters as enumerated below.   1. Number of features it's learning via auto-encoder   2. Weight decaying parameter of the auto-encoder   3. Weight decaying parameter of the logistic regression   4. Sparsity parameter of the auto-encoder   5. Weight of the sparsity penalty term of the auto-encode I’m planning to use random search for obtaining optimum values for these parameters. More information about random search for hyper-parameter optimization can be found in this paper My question is, in order to perform random search we need identify the appropriate ranges for each hyper-parameter. For example, the weight decaying parameter of the auto-encoder belongs to `[X,Y]`. So do you know a published paper to extract these initial values?


title: In R, find function F(x) to transform values in a vector to a normal distribution?

text: I have a PDF (Probability Density Function) generated from a vector of 1,000,000 empirical values. This empirical PDF is heavily skewed to the right. In this form, I can't make accurate predictions using a linear regression. To fix this, is there some method to find the function F(x) to transform (i.e. "squash") the values in the vector into a standard normal distribution, so I can feed said transformed vector into a linear regression? Of course, this would also involve finding the inverse of F(x) that transforms (i.e. "de-squashes") any predictions back into the original empirical PDF. **What I have tried** So far, I have managed to generate the density function from the empirical data: ![enter image description here](http://i.stack.imgur.com/HIBUP.png) Here is the R code:               par(mfrow=c(2,1))          install.packages("bootstrap")     library(bootstrap)     data(stamp)     nobs <- dim(stamp)[1]     hist(stamp$Thickness,col="grey",breaks=100,freq=F)      dens <- density(stamp$Thickness)     lines(dens,col="blue",lwd=3)          plot(density(stamp$Thickness),col="black",lwd=3, main="Simulation to choose density plot")      for(i in 1:10)      {       newThick <- rnorm(nobs,mean=stamp$Thickness,sd=dens$bw*1.5)       lines(density(newThick,bw=dens$bw),col="grey",lwd=3)     }          # If I wanted to do a linear regression to predict stamp thickness,     # what is the function F(x) to "squash" (i.e. transform) the "stamp"     # vector into a normal distribution, and the corresponding inverse      # function Finv(x) to "desquash" (i.e. untransform) any predictions back      # into the original prediction?      **Update 1** @Andre Silva sugggested that: > What need to have normal distribution are the residuals (predicted versus > observed) derived from your (multiple) linear regression model. According to post on Multiple Linear Regression: > After fitting the regression line, it is important to investigate the > residuals to determine whether or not they appear to fit the assumption of a > normal distribution. A normal quantile plot of the standardized residuals y > - is shown to the left. Despite two large values which may be outliers in > the data, the residuals do not seem to deviate from a random sample from a > normal distribution in any systematic manner. ![enter image description here](http://i.stack.imgur.com/3ybm0.gif) **Update 2** See Left skewed vs. symmetric distribution observed for R code that illustrates that the only relevant concern is if the residuals are normally distributed.


title: How to solve a dilemma case on Support Vector Machine?

text: I am learning the SVM classification and especially interested in applying to medical data. Now, I encounter a problem and do not know if this dilemma has a terminology for it. Assume that there are samples of healthy people ( of both gender) and people with liver cancer ( of both gender). If we label healthy people sample as class 1 and the people with cancer as class 2, we can train a binary SVM and obtain a classifier 1 to predict any new patient. Now, image another scenario. Assume that we first divide all samples by gender into two groups, one for female and the other for male. For each gender, we still label healthy patients vs cancerous patients into 2 classes and train a binary SVM to obtain classifier 2 and classifier 3 for female and male samples respectively. The question is if there is a new female patient, which classifier, 1 or 2, should be used to obtain more accurate prediction ? Here is the dilemma for the arguments I have (1) When the number of samples is large, the prediction should be more accurate. Based on this argument, the classifier 1 seems a good choice. (2) However, if we divide samples into female and male groups first, the classifier 2 seems a better choice since the new patient (unknown test sample) is female. Does this kind of dilemma have a terminology or does anyone know any further information or how to solve problem like this ? I am not even sure if this is a legit question and sorry for the naive question in advance. Thanks


title: Can I convert a statistical significance into a zscore using a normal distribution?

text: Suppose I want to turn statistical significance of 0.05 into a z-score of ~1.64 Can I use a normal distribution to convert these?


title: Models with Nested Effects

text: I am having difficulty understanding nested effects. Suppose 20 people each are assigned to treatment groups A B C and D. Would it be correct to say that Person is nested within treatment group? So for example, for each subject within each treatment group, each subject underwent 10 trials and for each trial, his heart rate was recorded. Would a model be: $HeartRate_{ijk} = Treatment_{i} + Subject_{j}(Treatment)_{i} + Trial_{k}(Subject_{j}) + e_{k(j(i))}$ where Subject and Trial are random effects and are nested? So there are no interaction terms in this model?


title: Phylogenetic dependent variables: ANOVA?

text: I understand deriving a covariance matrix from phylogenetic data to make $cov(X,Y) = 0$ for two variables you're making a regression on. But what happens if you have one continuous variable, that you've previously shown to be dependent on phylogeny, and one ordinal variable? The latter being ordinal, I'm not sure how to relate this to the way in which phylogenetic dependence results in biased test statistics. Is it meaningful to calculate Felsenstein's Phylogenetic Independent Contrasts on your continuous variable and use these for your ANOVA? The PIC value is: $$C_{ij} = \frac{(X_i - X_j)}{\sqrt{d_{ij}}} $$ Where $X_i$ is $X$ for species $i, X_j$ is $X$ for species $j$, and $d_{ij}$ is the pairwise distance between species $i$ and $j$ on the phylogenetic tree.