title: Wide-sense stationary linear process = ARMA process?

text: In Brockwell and Davis' Introduction to Time Series and Forecasting, a linear process is defined to be ![enter image description here](http://i.stack.imgur.com/FHA9V.png) An ARMA process/model is defined to be ![enter image description here](http://i.stack.imgur.com/7Hhg6.png) Note that by "stationary", the book means "wide-sense stationary". An ARMA process can also be written as a form of linear process ![enter image description here](http://i.stack.imgur.com/mtkKB.png)   1. Is it correct that a stochastic process is an ARMA process, if and only if it is a linear process and also wide-sense stationary?   2. Despite that the definition of an ARMA process says $X_t$ only depends on past values of $X$ and $Z$, is it correct that an ARMA process may not be causal, because the solution to (3.1.1) is (3.1.3)? Note a causal ARMA process is defined in (3.1.5) here. Do I misunderstand something? Thanks and regards!


title: What is the best one-step linear predictor for a white noise?

text: Suppose $X_i, i \in \mathbb N $ is a white noise. Given $X_1,\dots, X_n$, what is the best one-step linear predictor for $X_{n+1}$, best in the sense of mean square error? What if the question is changed to the following case: Suppose $X_i, i \in \mathbb N $ are uncorrelated, and they have identical mean and identical variance. Given $X_1,\dots, X_n$, what is the best one-step linear predictor for $X_{n+1}$, best in the sense of mean square error? Thank you very much!


title: How do I incorporate 2 peaks in a liquidity model using curve-fitting and nls

text: This is my first time at stack exchange, hence please pardon me if I miss something. I have multiple questions and I am tearing myself as curve-fitting and estimation is foreign land to me, I am ready to learn given directions and ideas. I am trying to model liquidity (ease of trading) of a security. In my case the security has very sparse trading. Once it is issued, trading is quite active and dies out then. After some years it may again have some hectic trading and becomes inactive for the rest of its life. Typically trading marketshare of the security depends on its age and outstanding amount. A sample datafile can be found here. In the data file ms_vol is the marketshare based on volume, ms_trds is the marketshare based on number of trades, osd is the outstanding amount and age is the age of the security. A pdf of ms_vol vs age plot can be found here. A research note has modeled such relationship as $$ MS = \beta_1 \exp[-\beta_2 (Age-\beta_3)^2] + \beta_4 \beta_5^{Age} $$ All the $\beta$'s are > 0. This models only one hump, not two. I tried to use nls in R unsuccessfully with the following code and error message               dfX = read.csv("trading.csv")     mod_msV = nls(       ms_vol ~ beta1*exp(-beta2*(age-beta3)^2) + beta4*beta5^age,       start=list(beta1=0.3, beta2=0.83, beta3=0.55, beta4=0.5, beta5=0.5),       data=dfX, trace=T)          29310.88 :  0.30 0.83 0.55 0.50 0.50      Error in numericDeriv(form[[3L]], names(ind), env) :        Missing value or an infinity produced when evaluating the model      Following are my questions:   1. How do I extend the above model to incorporate two humps which may occur some age apart?   2. The above model has only one independent variable, age. In my case osd is another variable I would like to incorporate. ms_vol behaves with osd just like age, though it may not show 2 humps.   3. How do I estimate teh parameters? The error message by nls baffles me completely. My hunch is it may be because of the starting values. I have tried several sets but all of them lead to same error message.   4. Somewhere I read about "Flowering Curve" and it may be a curve of choice. I have tried unsuccessfully locating info on it, but without success. Is there any other curve I can try out. Also how do I begin learning about this topic being a non-statistician? I am sorry for asking too many questions. I am at the edge of sanity. Please help me with this.


title: moderation analysis

text: I'm using demographic variables (age, number of children, gender, socioeconomic status, educational status) as moderating variables between Self differentiation (independent variable) and Marital Satisfaction (DV). I'm planning to use hierarchical regression is it correct? if in case I would not follow random sampling but instead use purposive sampling would that be possible? should I run tests for normal distribution? thank you very much for the response.


title: How to analyse separation of events in time?

text: I want to test whether the breeding periods of several closely related birds is significantly separated in time. What statistical test should I use? My data consists of several observations per species with the date of breeding initiation (i.e. egg laying dates). Could I just simply do an ANOVA on the Julian dates of breeding initiation? * * * I only have data of breeding birds, none of non-breeding ones. We only recorded breeding initiation of birds on nests. Each bird was recorded only once. I don't think survival analysis would fit this type of data. Am I wrong? (I'm not familiar with survival analysis). Breeding time of different species do not deviate from normality. Variances are equal among species. Here is what the data looks like. Each type and color of line represents a different species. ![Here is what the data looks like.](http://i.stack.imgur.com/GWDoj.png)


title: How to model influence of group on mean and variance of heart rate in patients with replicated measurements?

text: I'm trying to model correctly the data of a research were doing. We have many indicators measured along a day in n patients. For example, we have heart rate recorded once per minute for every patient. The patients belong to 4 groups. We want to estimate the effect of the group they belong to on the mean value and on the variability of the single patient heart rate. What would be a statistical correct way to approach this? Should i treat each patient as a point value of mean and variance of their data or should I use each measurement in a complex model? I would like to use a glm model, simple multivariate or more complex mixed effects, but i'm not sure how to handle them. Would it make sense to model the dependent variable as the single patient measurement in time or I should model directly the mean and variance for every patient? In the first case I should take in account the intra person variability, and separate it from it the effect of the group. From reading around I thought I should use a mixed effect model with the person within group as random effect and the group as fixed, but how to treat the fact that the measurement are linked IN the single person (we could say paired)? and then how can I analyze the variance in a person with this model? The second model is easier (one less dimension), since I would already have mean and variance calculated for every person, then I would just do a regression to asses how group influence mean and variance. Very simple but I would loose a lot of power (instead of hundreds of measurements per person I would just have 40 patients really unevenly distributed). What would you suggest? To be noted, I excluded ANOVA because the numerosity of the values for every patient varies and even more the number of patients per group. Secondly, every patient has a different distribution of values especially between groups, and this distribution could not be said normal (sometimes are even bi o multimodal) Ex, two subjects frm different groups, with values in time and distribution: ![enter image description here](http://i.stack.imgur.com/VbOzl.png) ![enter image description here](http://i.stack.imgur.com/77KEI.png)


title: What are the advantages of hidden Markov model over neural network?

text: I would like to know how about the differences between hidden markov model and neural network.


title: Transforming data to meet regression assumptions for SEM including interaction effects

text: If you are fitting a SEM and wish to add an interaction effect between two variables, I understand that you convert the variables to z-scores, multiply them together and insert this new variable into the model as the interaction effect. What do you do if the interaction effect variable does not meet assumptions for regressions (e.g. it is not normal, has significant kurtosis etc)—is this an issue? Can you attempt transformations of the interaction variable, or does it not have to meet the regression assumptions? If you did transform it, would this not create problems (i.e. you would be changing the interaction variable relative to the two variables that are already in the model)?


title: Better way to get best-fit regression formula (forcing y-int = 0)

text: I'm looking at a data set and trying to understand relationships between what is charged and what is paid. First pass with a simple regression using excel's data analysis regression tool gives the following, without forcing zero intercept: ![Residuals Plot](http://i.imgur.com/IzaJJyF.png) Second pass with forcing zero intercept: ![Residuals Plot](http://i.imgur.com/G8HG4Uq.png) I get better R squared values (0.82 vs 0.35) by forcing the zero intercept, but it looks like it is only because the residuals balance out. Is there a better way to run a regression here to get a best fit formula?


title: Determining conditional probability when there are two events but one condition and vice versa

text: “Determine probability that a purchase is paid for by debit card OR credit card given that the customer paid over $100 for it” [2 events and 1 condition] “Compute the probability that a child survives an accident given that the child is either 3 OR 4 years old.” [ 1 condition, 2 events] My question is, how is this handled? What formula should I use? What is the notation? The OR really throws me off because I know OR means addition rule and in conditional probability, one must use multiplication rule.


title: Difficulty in calculating joint probability

text: I have a table.                          CASH   CREDIT CARD   DEBIT CARD   TOTAL                           UNDER 20   0.09          0.03         0.04    0.16     20-100     0.05          0.21         0.18    0.44       OVER 100   0.03          0.23         0.14     0.4            TOTAL      0.17          0.47         0.36       1        The question is: “Find the probability that a debit card purchase was $20 and over”. In working this question, I became confused. I came up with several possible answers. I would like to know which one is right. How do I work out this question? Is the probability:   1. $0.18+0.14=0.32$   2. $ \frac{0.18}{0.44} + \frac{0.14}{0.40} = 0.41+0.35= 0.76$   3. $ \frac{0.18+0.14}{0.36} = 0.89$   4. $ \frac{0.14+0.18}{0.84}= 0.38$


title: Testing statistical significance of male vs female enrolment preference percentages between two different fields of education

text: The problem I have is as follows. I have enrolment data describing the number of male and female preferences for two broad fields of education: i) Information Technology and ii) Engineering. The data reflects "first preferences", that is, the number of people who wrote down on paper that their most preferred degree was either one in the field of Information Technology or the field of Engineering. The data does not reflect actual enrolments. The numbers are as follows: **Information Technology:**   Female: 266 (~13%)   Male: 1783 (~87%) **Engineering:**   Female: 684 (~12.5%)   Male: 4773 (~87.5%) The question I am trying to answer is: _is there a statistically significant difference between the gender ratios of Engineering versus Information Technology, and at what level of significance?_ I apologize if this is a noob question. I only have relatively basic background knowledge in statistics and tests of significance and this problem didn't seem to fit any of the methods I am familiar with. I'm hoping someone might be able to point me in the right direction as to which statistical test(s) would be appropriate given the dataset and any caveats I should be aware of.


title: Using doctor's data to identify hospitalisations

text: I have access to two large medical datasets of observational records in the UK. The first - Clinical Practice Research Datalink (CPRD) - has data on 100,000's of patients - largely dates of doctor's visits, diagnoses & prescriptions. Crucially, it has little information on hospitalisations. The second dataset is of hospitalisations - Hospital Episode Statistics (HES) - and importantly, I can identify HES patients in the CPRD. The HES data is a lot smaller - perhaps 5% of CPRD patients (I haven't the figures to hand)? - but it's still tens of thousands of patients. The brief is to look for signals in the CPRD data that might indicate a hospitalisation - using the HES data to verify this. I'll likely do this within a cohort of patients with a specific illness (respiratory disease, for example) which should make their medical patterns more homogeneous than the general population. The reason for the research is that CPRD is richer in medical history making it the most useful - but lacks that crucial element of hospitalisations. Still, it does have some - a patient might visit the doctor's and be told to go to ER, for example - this will be recorded in CPRD. Ultimately, I'm not optimistic I'll find much but I need to give it a good go. My background is traditional frequentist statistics where tools such as agreement measures just seem woefullly inadequate. I have a nodding knowledge of some bayesian and machine learning techniques - textbook rather than practical - and the latter seems appropriate in this case. I'm really looking for keywords & ideas to research (not necessarily full answers) as this will help greatly. Also, I'll likely program this in SAS 9.1.3 (or 9.3 eventually) - if that makes any difference - though I might be able to use R (though dataset size would definitely be an issue). I'm grateful for any help. EDIT: Resp Arthur Small: Yes, I suppose you could use CPRD records in that way, using HES records to 'confirm' hosps. I wonder though - given that CPRD is longitudinal data - at which points in CPRD you determine these probabilities (without first peeking at HES to get the date). Still, it definitely gives me something to investigate further - and is in my world of inferential stats! EDIT: Resp Peter Ellis: It is possible to use CPRD patients with HES (and the overlapping period) so yes, we can be sure that all hosps are recorded in HES for that subset of CPRD patients. I can't see how logistic regression would work however. If a patient has a HES hosp in June and Nov of a given year, I'm not sure how regressing on CPRD events will point to 1 hosp and not the other. I think it's the time element that's confusing me. I 'anticipate' a cluster of CPRD events (and event types) around HES hosps - it's that sort of analysis I believe I need. Thanks.


title: What does it mean when all edges in a real-world network/graph are statistically just as likely to happen by chance?

text: I've been using the backbone network extraction method outlined in this paper: http://www.pnas.org/content/106/16/6483.abstract Basically, the authors propose a method based in statistics that produces a probability, for each edge in the graph, that the edge could have happened just by chance. I use the typical statistical significance cutoff of 0.05. I've been applying this method to several real-world networks, and interestingly some networks end up with no edges as significant. I'm trying to understand what this entails for the network. The only other time I've applied the method to a network and had no edges come out as significant was when I applied the method to random networks that I generated, which is exactly what we'd expect. As an example real world network, you may have seen the recent network visualization that went on The Economist showing the polarization of the U.S. Senate in the past 25 years: http://www.economist.com/news/united- states/21591190-united-states-amoeba. I applied the backbone network extraction method to those networks and no edges came up as significant. Even though the raw edges apparently show preferential attachment and clustering, is this just by chance? Is the Senate voting network network essentially random?


title: How to decompose this conditional probability?

text: Could someone please tell me, where this comes from: $$ p(\beta, \sigma^2 | y, \tau) \propto p(y | \beta, \sigma^2) p(\beta | \tau) p(\sigma^2). $$ Thank you!


title: Testing which groups have had a significant change

text: I have three groups with two hundred individuals each, and each individual within a group has a pre and post treatment measurement. A single group contains individuals with similar pre measurement scores. If I wanted to see which groups have had a significant change in their measurements from pre to post treatment, what statistical test should I use?


title: p-values different for binomial test vs. proportions test. Which to report?

text: The following binomial test (exact binomial test) has a _p_ -value of **0.05698** , greater than 0.05 by a bit. 95% CI is (0.5698, 0.8077).               binom.test(x = 44, n = 63, p = 0.8, alternative = "two.sided")      The following proportions test (without Yates' continuity correction) has a _p_ -value of **0.04382** , less than 0.05 by a bit. 95% CI is (0.5764, 0.7976).               prop.test(x = x1, n = n1, p = 0.8, alternative = "two.sided", correct = FALSE)      So, which test should I choose? Which _p_ -value to report? Thanks.


title: Estimating the error in the average of correlated values

text: **tl;dr** I can only generate samples of a random variable $X$ using MCMC. How can I find the error in the estimate of the expected value of $X$ based on this MCMC data? * * * ### The problem I have a "black box" which gives me a sequence of "random" numbers (MCMC sampler). I have recorded $n$ numbers from this black box, which is all the data I can work with (I can't sample additional points). Now I need to estimate the average of the black box based on this data, and _get some idea about how precise the estimate is_. In other words, if I simply calculate the average of these $n$ datapoints, how far is it going to be from the true average of the black box? The difficulty: Consecutive numbers coming form the black box are strongly correlated, however the $i^\text{th}$ and $(i+k)^\text{th}$ numbers are effectively uncorrelated if $k$ is large. Are there any methods to deal with this type of situation, and find the error in the estimate? Let's assume that $n$ is large enough that I have some uncorreltaed data points. * * * ### What I'm doing now In the following I'll describe what I am doing right now and why I am not satisfied with this method. This method works but it is not ideal. Let's look at the case when the numbers $x$ are uncorrelated. Let $\langle x \rangle_n$ denote the average of $n$ consecutive numbers, which is an estimate of the expected value of the black box. Then $\operatorname{Var}(\langle x \rangle_n) = \operatorname{Var}(x)/n$. $\sqrt{\operatorname{Var}(\langle x \rangle_n)}$ can be used as the error of the estimate $\langle x \rangle_n)$ and it's easy to compute based on $\operatorname{Var}(x)$. However my actual data is correlated. So I partitioned my dataset of $n$ elements into runs of $k$ elements and calculated the product $k \times \operatorname{Var}(\langle x \rangle_k)$ for various $k$. When $k$ is large enough, and the correlations are not important any more, this product will settle at a constant value, which can be used to extrapolate for $\operatorname{Var}(\langle x \rangle_n)$. This is a plot of $k \operatorname{Var}(\langle x \rangle_k)$ for various values of $k$ for four different datasets (four colours) of size $n = 2^{27}$: ![enter image description here](http://i.stack.imgur.com/k4rRY.png) You can see that as $k$ becomes large enough, the product $k \operatorname{Var}(\langle x \rangle_k)$ settles on a constant. But as $k$ becomes too large, I can only parition my length-$n$ dataset into a few parts of length $k$, so the estimate of the variance becomes inaccurate. Hence the noise on the right of the curve. So I need to find the flat part of the curve and discard the noisy part, then estimate the variance of the average of the full dataset ($\langle x \rangle_n$) based on it. ### Why my method isn't ideal The method I describe above is complicated and involves several steps, and several choices in each step. For how many (and what magnitude of) $k$ values should I compute the variances? How do I determine where is the flat part of the curve? How do I decide which part is too noisy to keep? All of these are separate problems on their own with many valid heuristic solutions, which will all influence the result. Also, this method isn't very efficient computationally. I am hoping for a simpler and cleaner robust method where the details won't have such a large effect on the end result.


title: Why did removing the between-subjects factor from a two-way mixed ANOVA change the SSwithin?

text: I have been playing with this SPSS data file, which is a class example from a university I am not enrolled in. In the original file there are 10 men and 10 women. I did a two-way mixed ANOVA with gender as the between-subjects factor and I also used a within- subjects factor "looks". Then I got rid of gender and just used the within- subjects factor in a within-subjects ANOVA. The SS for looks was the same either way. Then I deleted two women from the sample and undertook the same analyses again. Now the SS for looks in is no longer the same in the two analyses. Why does this happen? ![enter image description here](http://i.imgur.com/u7aViTo.png)


title: Fixed effect model is not significant without year dummies

text: My fixed effect model is insignificant when I don't include the year dummies $$\text{Prob $>$ F} = 0.7769$$ When I try other regressions the model is fine, for example Driscoll & Kraay fixed effect. When I delete some outliers from my dependent variable, the fixed effect become significant and the random effect model turned insignificant! Is there a problem in the model of any variable? What's your opinion?


title: K-fold Cross Validation and Training/CV/Test set Techniques for choosing regularization parameter of Regression

text: Suppose I want to fit a lasso/ridge regression to a training set. Then, I need to choose $\lambda$, the regularization parameter. To choose $\lambda$, I can use two methods:   1. K-fold Cross Validation (from An Introduction to Statistical Learning p. 227):     1. Divide the training set into K folds (randomly).     2. Choose one fold, fit with the data in K-1 folds. Do it for all K folds.     3. Average the error for each model.     4. Choose $\lambda$ that gives the lowest error.     5. Fit with all the data in original training set with $\lambda$ chosen in (4).   2. Training/Cross Validation/Test Sets method (as taught by Andrew Ng in Coursera):     1. Divide the original training set (randomly) into 3 subsets, (new) training set, cross validation set, and test set, with proportion approx. 60%, 20%, 20%.     2. Fit with the new training set for every value of $\lambda$ you determined.     3. Measure the error of the model using cross validation set for each $\lambda$.     4. Choose the model and $\lambda$ which gives the lowest cv error.     5. Test the model with $\lambda$ chosen at (4) on test set to measure the error. Which of these methods gives lowest bias and variance when fitting lasso/ridge regression?


title: Hypothesis testing question, references

text: I'd appreciate reference books/online resources in learning how to apply A/B hypothesis testing for problems of the following type: Given an initial strategy B with $X$ samples yields a success rate of $b$. Now say we are given two new strategies $B_1$ and $B_2$, one trial each with number of samples $X_1$ and $X_2$, and success rates $b_1$ and $b_2$, respectively. Is strategy $B_1$ better than $B_2$, with respect to the success rate $b$ in strategy B? Further, in using measures like $p$-values, how much does our decision of $B_1$ vs $B_2$ change if we have now, say 1000 trials of both $B_1$ and $B_2$ each (so now a total of $1000\cdot X_1$ samples for $B_1$ and $1000\cdot X_2$ for $B_2$)? Aside: how might the Neymann-Pearson lemma apply and help us with this type of decision question _in practice_?


title: Statistically significant with binomial

text: What does this statement mean: "we calculated binomial 95% confidence intervals to determine the number of weeks that were statistically- significantly outside of the bounds of error"? How can one calculate it?


title: Choice of plot or statistic for data with wide range

text: I have some simple datasets that consist of a discrete integer-valued variable together with frequency counts. The frequencies vary a lot as does the value of the variable. For example one set has values from 1 to a little less than 15000 with frequencies from 1 to half a million. The small values have the many occurrences, e.g. no value higher than 10 occurs more than 50000 times whereas the values 1 through 10 together appear around 1.5 million times. I am looking for a good kind of plot (preferably) or statistic to use for presenting and overviewing this data. What is a good choice? At the moment I am plotting a line segment from (n,0) to (n, freq) for the data point (n, freq) but due to the distribution this is not very informative.


title: Error with RLRsim package while testing for variance components

text: I am trying to fit the mixed effects model and test for a variance component using RLRsim package. $$ Y_{ik} = \alpha + \mu_{k} + \beta G_{i} + Z_{k}G_{i} + b_{i} + \epsilon_{ik} $$ Random effects: $$ Z_{k} \sim N\left(0,\psi\right) \qquad b_{i} \sim N\left(0,\tau\right) \qquad \epsilon_{ik} \sim \left(0,\lambda\right) $$ I am testing for the variance component of the random effect $Z_{k}$ in a constrained setting, where $H_0 : \psi = 0$ and $H_A:\psi > 0$ using RLRsim package. Here's how my data looks like --                   > head(data,14)           IND   Gene Tissue       Geno       1    1  5.503      1 -0.8952456       2    1  7.838      2 -0.8952456       3    1  6.903      3 -0.8952456       4    1  9.821      4 -0.8952456       5    2 12.276      1  0.6482813       6    2 16.922      2  0.6482813       7    2 15.709      3  0.6482813       8    2 17.116      4  0.6482813      ...with a total of 150 observations and 4 tissue types.                   m = lmer(Gene ~ Geno + Tissue + (0+Geno|Tissue),data,REML=T)         m0 =lmer(Gene ~ Geno + Tissue + (1|IND),data,REML=T)         mA =lmer(Gene ~ Geno + Tissue + (1|IND) + (0+Geno|Tissue),data,REML=T)         rlrt = exactRLRT(m=m, mA=mA, m0=m0, nsim=10000)      However, I am getting the following error message when I run exactRLRT().                   Error in chol.default(cov2cor(Vr)) : the leading minor of order 2 is not positive definite      Here is the summary of the model with the random effect that is needed to be tested --                   > summary(m)           Linear mixed model fit by REML            Formula: Gene ~ Geno + Tissue + (0 + Geno | Tissue)            Data: data            AIC  BIC logLik deviance REMLdev           2105 2136  -1045     2077    2091           Random effects:           Groups   Name Variance Std.Dev.           Tissue   Geno 0.000    0.0000             Residual      1.881    1.3715             Number of obs: 600, groups: Tissue, 4                Fixed effects:                        Estimate Std. Error t value           (Intercept)  7.84669    0.11198   70.07           Geno         5.00645    0.05618   89.12           Tissue2      4.23212    0.15837   26.72           Tissue3      4.10939    0.15837   25.95           Tissue4      4.34121    0.15837   27.41                Correlation of Fixed Effects:                      (Intr) Geno   Tissu2 Tissu3             Geno     0.000                                  Tissue2 -0.707  0.000                           Tissue3 -0.707  0.000  0.500                    Tissue4 -0.707  0.000  0.500  0.500      Here's the summary of the model under the null --                   > summary(m0)           Linear mixed model fit by REML            Formula: Gene ~ Geno + Tissue + (1 | IND)            Data: data            AIC  BIC logLik deviance REMLdev           1940 1970 -962.8     1911    1926           Random effects:           Groups   Name        Variance Std.Dev.           IND      (Intercept) 0.92302  0.96074            Residual             0.96262  0.98113            Number of obs: 600, groups: IND, 150                Fixed effects:                        Estimate Std. Error t value           (Intercept)  7.84669    0.11212   69.99           Geno         5.00645    0.08837   56.65           Tissue2      4.23212    0.11329   37.36           Tissue3      4.10939    0.11329   36.27           Tissue4      4.34121    0.11329   38.32                Correlation of Fixed Effects:                   (Intr) Geno   Tissu2 Tissu3          Geno     0.000                               Tissue2 -0.505  0.000                        Tissue3 -0.505  0.000  0.500                 Tissue4 -0.505  0.000  0.500  0.500      Here's the summary of the alternative model --                   > summary(mA)           Linear mixed model fit by REML            Formula: Gene ~ Geno + Tissue + (1 | IND) + (0 + Geno | Tissue)            Data: data            AIC  BIC logLik deviance REMLdev          1942 1977 -962.8     1911    1926          Random effects:          Groups   Name        Variance Std.Dev.          IND      (Intercept) 0.92305  0.96075           Tissue   Geno        0.00000  0.00000           Residual             0.96262  0.98113           Number of obs: 600, groups: IND, 150; Tissue, 4               Fixed effects:                       Estimate Std. Error t value          (Intercept)  7.84669    0.11212   69.98          Geno         5.00645    0.08837   56.65          Tissue2      4.23212    0.11329   37.36          Tissue3      4.10939    0.11329   36.27          Tissue4      4.34121    0.11329   38.32               Correlation of Fixed Effects:                  (Intr) Geno   Tissu2 Tissu3          Geno     0.000                               Tissue2 -0.505  0.000                        Tissue3 -0.505  0.000  0.500                 Tissue4 -0.505  0.000  0.500  0.500      Here's the session info--                  > sessionInfo()          R version 3.0.1 (2013-05-16)          Platform: x86_64-apple-darwin10.8.0 (64-bit)               locale:          [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8               attached base packages:          [1] stats     graphics  grDevices utils     datasets  methods   base                    other attached packages:          [1] BiocInstaller_1.10.2 RLRsim_2.1-0         mgcv_1.7-24                   [4] car_2.0-18           nnet_7.3-7           MASS_7.3-27                   [7] lme4_0.999999-2      Matrix_1.0-12        lattice_0.20-15                    loaded via a namespace (and not attached):          [1] cluster_1.14.4    gdata_2.13.2      gplots_2.11.3     grid_3.0.1                 [5] gtools_2.7.1      Hmisc_3.10-1.1    lmerTest_1.2-1    nlme_3.1-109               [9] numDeriv_2012.9-1 pbkrtest_0.3-5    stats4_3.0.1      tcltk_3.0.1               [13] tools_3.0.1            I appreciate any help here. Thanks!


title: Non-informative prior for integer restricted parameters

text: Does anyone have advice/references on choosing a minimally informative prior for the shape parameter of an Erlang distribution (gamma distribution with shape parameter restricted to integer space), or a joint prior on both parameters.   My gut tells me that a discrete uniform could be a bad idea. A truncated poisson would work but will be quite informative....


title: Calculating proportions by age in R

text: I working on this data set with marital status and age. I want to plot the percentage of never married man versus each age. Could you please help me to figure out the way how to do it in R? So far I have created two separate arrays with men never marries and ever married. I know how many case of each I have. What i need to do is to count number of people that were never married at each age and divide it by the total number of never married people to get a percentage. I hope I was clear. Thank you


title: The homogeneity of error variance in MMR with continuous moderator?

text: I was wondering if I have to worry about the homogeneity of error variance assumption (which states that the predicted scores for Y should be similarly distributed about the regression line for each of the moderator based populations) when my moderator is continuous? I know you should not dichotomize a continuous moderator but that is the only way I would be able to get "moderator-based" subgroups. I have found literature on how to test this when the moderator is categorical and I am not sure if it is something I have to worry about. If I do have to check this assumption how would I go about it? Thanks.


title: How to fit a multiple regression model that takes into account 50 days of each individual's history

text: I am using simple linear regression to predict the number of days of school a single student might miss after 250 days using their cumulative daily attendance data over the first 50 days of the student's schooling--i.e., predicting how many days of school they will have missed at 250 days based on their individual cumulative 50 day trajectory. This works fairly well for about 75% of students. As a next step, I want to refine Y = a + bX to test whether I can more accurately predict the number of days a student might miss after 250 days (based on their first 50 days), taking into account their gender, ethnicity, etc., using multiple regression. I'd understand how to do this if I were using regression to predict an outcome based on a group of single co-ordinates, i.e. predicting a test result based on the number of days missed for 100 males / females. But, I'm not sure how to generate a regression model that takes into account each 50 day matrix of co- ordinates for all 100 student simultaneously. I then want to take the general effects of gender, ethnicity on the slope of this model (i.e., Y = a + b1X1 +b2X2...) and plug it back into the individual student's regression line to test if I can more accurately predict their 250 days outcome. Could anyone offer me some starting guidance on conducting regression in this situation please? Or indeed even if regression is the appropriate choice? I'm using Oracle Analytics for student linear regression and was planning to use SPSS for multiple regression.


title: Why is the Bonferroni procedure applicable only when the effects to be investigated are identified in advance of the study?

text: From Kutner's Applied Linear Statistical Models > In the next three sections, we shall discuss three multiple comparison > procedures for analysis of variance models that permit the family confidence > coefficient and the family a risk to be controlled. > > Two of these procedures, the Tukey and Scheffe procedures, allow data > snooping to be undertaken naturally without affecting the confidence > coefficient or significance level. > > The other procedure, the Bonferroni procedure, is applicable **only when the > effects to be investigated are identified in advance of the study**. Why is the Bonferroni procedure "applicable only when the effects to be investigated are identified in advance of the study"? Is it the same case for the Tukey and Scheffe procedures? Thanks!


title: For the model $Y = α + βx + \epsilon$, estimate mean and predicted value of Y when $x = x_0$

text: Consider n independent observations {($x_i, y_i$) : $1 \le i \le n$} from the model $$Y = α + βx + \epsilon,$$ where $\epsilon$ is normal with mean 0 and variance $σ^2$. Let $\hat α$, $\hat β$ and $\hat σ^2$ be the maximum likelihood estimators of $α$, $β$ and $σ^2$ respectively. Let $v_{11}$, $v_{22}$ and $v_{12}$ be the estimated values of Var($\hat α$), Var($\hat β$) and Cov($\hat α$,$\hat β$) respectively. Then (a) What is the estimated mean of $Y$ when $x = x_0$? Estimate the mean squared error of this estimator. (b) What is the predicted value of $Y$ when $x = x_0$? Estimate the mean squared error of this predictor.


title: How to statistically analyse the survey data of the same sample taken at two points in time

text: I am conducting a study on a sample of 100 people in order to test whether or not education has an effect on employment status. Basically i'm running a logistic regression where education is the independent categorical variable (i.e 1= HighSchool, 2=Bachelors, 3=Masters etc...) and the dependent variable being binary response (1= employed, 0=unemployed). I also conduct this survey again 4 years later with the exact same group of people. Hence the variable that has changed is age (meaning people in my sample are 4 years older). How can I analyse this type of study? what are some of the tests I can conduct? Please be as detailed as possible, i do not have a strong statistics background and I would greatly appreciate any guidance you can provide :) Thanks.


title: Diagnostics for GEE in R

text: I have been checking out which diagnostics to use for a GEE analysis. It seem that influence measures are appropriate (Preisser, 1996). Does anyone know of a package that can be used in R to examine influence measures in GEE? The "stats" package does not work: it returns an error when I run               influence.measures(model)       on a gee model, but works fine with an ordinary linear regression.


title: Derive P(C | A+B) from Cox's two rules

text: I am working my way (self-study) through E.T. Jaynes' book _Probability Theory - The Logic of Science_ ## Original Problem _Exercise 2.1_ says: "Is it possible to find a general formula for $p(C|A+B)$ analogous to [the formula $p(A+B|C)=p(A|C)+p(B|C)-p(AB|C)$] from the product and sum rules. If so, derive it; if not, explain why this cannot be done." ### Givens The rules I have to work with are: $p(AB | C) = p(A|C)p(B|AC) = p(B|C)p(A|BC)$ and $p(A|B)+p(\bar{A}|B)=1$ Where we can also use logical identities to manipulate propositions. For example: $A+B=\overline{\bar{A}\bar{B}}$ ### Assumption of Solvability I believe it must be possible because he does not introduce any other rules later and having a simple logical combination of propositions that was not easily expressible would defeat Jaynes' central thesis. However, I've been unable to derive the rule. ## My Attempt To keep myself from getting confused due to using the same variable names as the givens, I am solving the problem as: Derive a formula for $p(X|Y+Z)$ ### Introducing a tautology for conditioning My best attempt at solving it so far has been to introduce a proposition $W$ which is always true. Thus I can rewrite $Y+Z$ as $(Y+Z)W$ (since truth is the multiplicative identity). Then, I can write: $p(X|Y+Z)=p(X|(Y+Z)W)$ So, rewriting one of the givens as Bayes' rule: $p(A|BC)=\frac{p(B|AC)p(A|C)}{p(B|C)}$, I can write: $p(X|(Y+Z)W)=\frac{p(Y+Z|XW)p(X|W)}{p(Y+Z|W)}=\frac{p(Y+Z|X)p(X|W)}{p(Y+Z|W)}$ ### Why this doesn't work The term $p(Y+Z|X)$ is easy to deal with. (Its expansion is referred to in the problem definition.) However, I don't know what to do with $p(X|W)$ and $p(Y+Z|W)$. There is no logical transformation I can apply to get rid of the $W$, nor can I think of any way of applying the given rules to get there. ## Other places I've looked I've done a Google search, which turned up this forum page. But the author does the same thing I tried without seeing the difficulty I have with the resulting conditioning on the introduced tautology. I also searched stats.stackexchange.com for "Jaynes" and also for "Exercise 2.1" without finding any useful results.


title: Specifying selectively-nested random factors in R lmer

text: Partly as a learning exercise, I’m attempting to reproduce a published SAS PROC MIXED analysis using R’s lme4 package. A full description of both the application and the statistical method is available in http://www.biomedcentral.com/1472-6750/11/15 . Briefly, the broad purpose is to assess the food safety of a specific genetically-modified crop (GMO). In the experiment, the GMO crop is grown along with a genetically similar but unmodified comparator and with a diverse set of commercially-available reference varieties. The reference varieties are presumed safe and thus provide estimates of the range of acceptable nutrient levels. The crop is harvested and subjected to chemical analyses to measure various constituents. Nutrient levels are then compared between the GMO and the comparator and between the GMO and the set of reference varieties. The published method defines a variety group variable with 3 levels: GMO, comparator (comp), and reference (ref). An integer index is also defined that singles out the reference varieties. To clarify, assume that we have 4 reference varieties. The variety-related variables are then:               variety vgroup  isref     comp    comp    0     GMO     GMO     0     refA    ref     1     refB    ref     1     refC    ref     1     refD    ref     1      Focusing only on this part of the statistical model and ignoring the other (important but straightforward) model factors (i.e. locations & blocks), the published SAS PROC MIXED formulation is:                   CLASS variety vgroup …;         MODEL y=vgroup …;         RANDOM isref*variety …;      Modeling vgroup as a fixed factor enables formal tests of (GMO – comp) and (GMO – ref). Note that isref is not designated as a CLASS variable and thus remains a numeric quantity. The SAS model runs successfully. A single variance is estimated for the isref*variety term, indicating that this quantity is a single factor rather than an interaction. The random effect predictions for the isref*variety factor are identically zero for the comp and GMO varieties and are non-zero for the ref varieties. I have yet to replicate these results in R lmer. The following model statement                   y ~ vgroup + (0 + factor(isref)|variety) + …      is close, but doesn’t annihilate the isref == 0 variance estimate and the associated random effect predictions. Perhaps as a direct consequence, the vgroup standard errors for the comp and GMO levels are considerably larger than the corresponding SAS results. So finally, is there a way to specify this exact model in R lmer? If so, how? If not, why?


title: Means of groups A and B differ significantly. I want to classify values into A or B

text: My data look like this. ![Heights of hoverers and non-hoverers](http://i.stack.imgur.com/Qi08m.png) The variable on the $x$ axis is height in inches. The variable on the $y$ axis is whether someone hovers when urinating at a public toilet. Each of the 103 points is a female university student. According to a $t$-test, the mean height of hoverers is significantly higher than that of non-hoverers. I would like to make an inference in the opposite direction; I would like to know whether taller women are more likely to hover.   1. Is it reasonable to suspect based on my $t$-test that the opposite inference would hold?   2. A logistic regression found no significantly increased odds of hovering corresponding to an increase in height. Is there a different test I should try for this opposite relationship?


title: Is there a known MLE for the numerator df of a sample of F statistics?

text: Suppose you have observations $f_1, f_2, \ldots, f_n$ which are drawn _i.i.d._ from a central F distribution, with unknown numerator degrees of freedom, $n_1$ and known denominator degrees of freedom, $n_2$. What is the preferred method for computing/estimating the MLE of $n_1$? If the likelihood is known to be concave as a function of $n_1$ (I would guess it is, but am not certain), I could just use a numerical method like golden section search. Is there anything better? I am aware of Spruill's method for ML estimation of the non-centrality parameter when the $f_i$ are noncentral, and $n_1$ and $n_2$ are known. It appears to use the transformation $$w_i = \frac{f_i}{n_2/n_1 + f_i}$$ Is this a standard transformation for F distributed RVs that might be applicable to this problem, or just a red herring?


title: Using hypothesis-testing to identify misleading picture captions

text: I have developed a web application in which I show a picture and two captions. One of those captions is the correct and the other is randomly selected from a corpus of picture captions. Next, I present the following question: "Which caption best describes the picture?" and the user casts a vote on the caption that he chooses. Since I store the user votes, I know how many votes each picture received and also, the fraction of votes that were on the correct caption. Since the 'correct' captions of some pictures are misleading, it is expected that, in these cases, the correct vote fraction is near %50 (that is, it is impossible to identify the correct caption). However, some pictures will have very descriptive captions, and the correct vote fraction is expected to be somewhat higher. Is it possible, given a confidence value (e.g. 95%), to identify pictures whose correct caption is impossible for the users to identify (correct vote fraction = %50)?


title: A/B/C testing evaluation over time ans across groups

text: I'm needing to validate an inherited testing project correct use of confidence testing. Direct mail creatives are being evaluated at month end over a 90 day campaign cycle for several unique customer groups (eg.: new, upgraded, lapsed, etc.) Here's a low volume data example for just one customer group (new) for a single month (Note, the eval process repeats for the next 2 months as the #mail & %Resp(onse) accumulates):               Treatment       #Mail #Response %Resp  Z-score       ---------       ----- --------- ------  -------       Control          182    35      19.23%   N/A       Treatment A      180    45      25.00%  1.33       Treatment B      189    28      14.81% -1.13       Treatment C      188    61      32.45%  2.94      The question (s):   1. Is this the most appropriate (single tail) test given required assumptions for normality & randomness of the data?   2. 'Treatment C' w/ a Z-score of 2.94 suggests a 95% confident lift over Control for month 1, declines to non-significance in month 2, but then reverts back to 95% confidence level in month 3. Should 'Treatment C' be considered the NEW control for the upcoming next quarter?   3. Now assuming a dissimilar performance for 'Treatment C' is also observed across the other remaining customer groups (upgraded, lapsed, etc.) during the same monthly eval process over the quarter, how should the analyst use that finding comparatively in determining whether 'Treatment C' should be used at all for the new customer group?   4. Is there a less error prone/subjective approach to solving this problem of what alternative creative should become the new control within and across groups over time?   5. Based on #4 above, is their a simple library or query example one can implement this within MySQL versus currently within Excel?


title: How to get started in data-relation algorithms and mathematics?

text: I am very interested in the concepts and discussions taking place here. However, I'm not entirely sure what this field of study is called or what the many branches of study are being discussed here are. It all seems to be related to finding relations and patterns in data so we can use it in a meaningful way. As a developer, I find myself increasingly researching things like Statistical Graphing, Natural Language Processing, Bayesian filters or Probably, Bin- Packing problems, and many more. Unfortunately, I'm not sure where to go to learn these things and better contribute to the community here. What names can I use to research and learn the topics discussed here? What would the equivalent college courses for this type of study be called?


title: How to detect a number of binomial distributions from a set of data?

text: Say I have n sets of data where each set of data has a binomial distribution. The n sets of data are then combined into one set. If I take the combined set of data, how would I go about detecting that this set is made up from n individual sets ? I have only a limited background in statistics, so here is an example to clarify the question: I collect the amount of energy used by an appliance at 3pm each day for a number of months Generally the amount of energy used follows these rules: - On Monday to Friday, the appliance uses 3000 kWh (mean value) - On Saturday, the appliance uses 200 kWh (mean value) - On Sunday, the appliance uses 50 kWh (mean value) I want to be able to take all the 3pm readings every day for a year and then be able to detect that the readings can be broken into the 3 sets as described above. After I deduce that there are 3 sets of data, I plan to calculate the mean value with std dev etc. for each individual set.


title: Excel's CHIDIST function in MATLAB

text: How to reproduce MS Excel's CHIDIST function in MATLAB? > Returns the one-tailed probability of the chi-squared distribution. The χ2 > distribution is associated with a χ2 test. Use the χ2 test to compare > observed and expected values. For example, a genetic experiment might > hypothesize that the next generation of plants will exhibit a certain set of > colors. By comparing the observed results with the expected ones, you can > decide whether your original hypothesis is valid. > > Syntax > > CHIDIST(x,degrees_freedom) > > X is the value at which you want to evaluate the distribution. > > Degrees_freedom is the number of degrees of freedom. So, if I have column of numbers as `x` and use the formula like `=CHIDIST(A1:A60,3)`, it returns 1 number. I have tried to follow two advices from this thread, but I always get a vector of numbers, not one number. Thanks


title: Estimating the distribution of a very large population of known size and unknown variance

text: I would like to estimate the distribution of a very large population of known size but unknown mean and variance. I cannot assume anything about the shape of the underlying distribution (although I am relatively certain that it is not normal). However, I am certain that the values in the population are non- negative, non-zero integers. I believe that the distribution is naturally lumped in some way discretely rather than being continuously distributed over the entire range. I cannot sample the entire population but I would like to estimate the probability density function. I would also like some level of assurance of the correctness of the estimated distribution. What is the best way to go about this?


title: Why multiplicative Holt-Winters requires strictly positive data points?

text: I've seen that multiplicative Holt-Winters requires strictly positive data points. I was wondering why it does not allow zero values?


title: Determining best fitting curve fitting function out of linear, exponential, and logarithmic functions

text: ### Context: From a question on Mathematics Stack Exchange (Can I build a program), someone has a set of $x-y$ points, and wants to fit a curve to it, linear, exponential or logarithmic. The usual method is to start by choosing one of these (which specifies the model), and then do the statistical calculations. But what is really wanted is to find the 'best' curve out of linear, exponential or logarithmic. Ostensibly, one could try all three, and choose the best fitted curve of the three according to the best correlation coefficient. But somehow I'm feeling this is not quite kosher. The generally accepted method is to pick your model first, one of those three (or some other link function), then from the data calculate the coefficients. And post facto picking the best of all is cherry picking. But to me whether you're determining a function or coefficients from the data it is still the same thing, your procedure is discovering the best...thing (let's say that which function is -also- another coefficient o be discovered). ### Questions:   * Is it appropriate to choose the best fitting model out of linear, exponential, and logarithmic models, based on a comparison of fit statistics?   * If so, what is the most appropriate way to do this?   * If regression helps find parameters (coefficients) in a function, why can't there be a discrete parameter to choose which of three curve families the best would come from?


title: Calculation of natural cubic splines in R

text: I am new to the use of cubic splines for regression purposes and wanted to find out 1) What is a good source (besides ESL which I read but am still uncertain) to learn about splines for regression?   2) How would you calculate the basis of a given natural cubic spline solution on new data? Specifically if one were to do the following:               data(iris)     colnames(iris)     Sepal.Length.ns<-ns(iris$Sepal.Length,df=5)     Sepal.Length.ns      How would you take the information in Sepal.Length.ns (knots, boundaries) and compute the values for a new observation? The reason is to code this process outside of R, once fit in R initially (i.e. to put a regression model using cubic splines into a production system). For example I can do this in R, but want to understand the calculation:               #three new observations to predict     newVector<-c(4.45,3.35,2.2)     pred.new<-predict(Sepal.Length.ns,newVector)      Thanks!


title: Binomial data analysis with all 0 responses for some treatment groups

text: We have a series of experiments where we measure virus transmission to plants when exposed to virus-infected insects for different time periods, so all of the experiments have similar types of independent and dependent variables. In one experiment, there are 6 time periods (1 to 24 hours) and 25 plants were tested (individually) for each time period. The response for each plant is yes or no (Plants are scored as virus infected). For 2 of the time intervals, all of the plants were negative for virus infection (0/25 for each time interval). I am using PROC GLIMMIX in SAS for the analyses. For all of the other experiments, using a binary distribution in the model statement gives reasonable results. For the experiment where two of the time intervals had 0 positive plants, if I use a binary distribution in the model statement, the standard errors for the two groups with 0 transmissions are huge, thus distorting the results. If I use a negative binomial distribution (based on counts of virus positive plants) the results seem reasonable. Since the same number of plants were tested for each time interval, using this approach works, but it differs from the other experiments. Is there a method to adjust/account for the zeros in treatment groups that would allow the binary distribution return reasonable results?


title: Determining the relationship between salesperson and products sold?

text: ### Context: I have a series of figures for car sales that show me (a) the usual number of car sales for particular models and (b) the number of car sales by a particular car salesperson for each model. Let's say the values are:               Model | Sales by all Salespeople | Sales by Salesperson x         A |                      100 |                     20         B |                       50 |                     40         C |                       50 |                      0      I want to find out if Salesperson x is significantly more responsible for sales of a particular model. My first thought is to determine the rate of sales per model vs the rate of sales per model for the salesperson, i.e.:               Model | % Total Sales for all Salespeople | % Total Sales for Salesperson x         A |                               50% |                             33%         B |                               25% |                             66%         C |                               25% |                              0%      However, naturally every salesperson has a significant variation in the cars they sell. ### Question:   * How do I determine if salesperson variation from the mean is statistically significant? ### Initial Thoughts: I _think_ the Pearson correlation has some bearing, as perhaps may chi-square distribution, but I really don't have the background yet to understand why, so introductory help is appreciated.


title: Probability for finding a double-as-likely event

text: Repeating an experiment with $n$ possible outcomes $t$ times independently, where all but one outcomes have probability $\frac{1}{n+1}$ and the other outcome has the double probability $\frac{2}{n+1}$, is there a good approximate formula for the probability that the outcome with the higher probability happens more often than any other one? For me, $n$ is typically some hundreds, and $t$ is chosen depending on $n$ such that the probability that the most likely outcome occurs most often is between 10% and 99.999%. In the moment I use a small program that calculates a crude approximation by assuming that the counts for how often each outcome shows up in $t$ trials are independent and approximate the counts using the Poisson distribution. How can I improve on this? **EDIT:** I'd strongly appreciate comments/votes on the two (maybe soon more) answers given. **EDIT 2:** As none of the two answers is convincing me, but as I don't want to let the 100 points bounty to vanish (and as nobody voted for/against one of the two answers), I'll just pick one of the answers. I'd still appreciate other answers.


title: Differences-in-Differences coefficients meaning (from Mostly Harmless Econometics)

text: In _Mostly Harmless Econometrics_ , section 5.2.1 (Regression DD), pages 233-234, equation (5.2.3) defines $Y_{ist}=\alpha + \gamma NJ_{s}+\lambda d_{t}+\delta (NJ_{s}.d_{t})+\epsilon_{ist}$, where $NJ_{s}$ is a dummy denoting observations from New Jersey and $d_{t}$ is a dummy for observations obtained in November. Considering $E[\epsilon_{ist} \vert s,t]=0$, the books states that: $\alpha=E[Y_{ist} \vert s=PA, t=Feb] = \gamma_{PA} + \lambda_{Feb}$ $\gamma=E[Y_{ist} \vert s=NJ, t=Feb] - E[Y_{ist} \vert s=PA, t=Feb] = \gamma_{NJ} - \gamma_{PA}$ $(...)$ Although, in **my opinion** : $E[Y_{ist} \vert s=NJ, t=Feb]=\alpha+\gamma$, so that $E[Y_{ist} \vert s=NJ, t=Feb] - E[Y_{ist} \vert s=PA, t=Feb] =\alpha+\gamma -\alpha=$ $=(\gamma_{PA}+\lambda_{Feb}+\gamma_{NJ})-(\gamma_{PA}+\lambda_{Feb})=\gamma_{NJ}$, where $\gamma_{NJ} \equiv \gamma$ Which has a different meaning from what the book stated (cited above). Could anyone explain what I am missing or if there is a mistake in the book? Found nothing about it at the blog.


title: plot for a binomial glm

text: Can any one help? I am trying to graph the probability of predicting a male bird based on tarsus length according to the GLM. I have run a GLM and I am not sure why I can't plot the predictions. My code and data follows.               Tarsus   Sex     33.49   female     36.36   female     36.61   female     37.27   female     37.95   female     38.39   female     38.5    female     39.5    female     40.11   female     40.3    female     40.39   female     40.47   female     41.14   female     41.4    female     41.5    female     41.94   female     43.64   female     41.4    male     41.55   male     41.8    male     41.98   male     42.51   male     42.56   male     42.7    male     43.04   male     43.1    male     43.39   male     43.61   male     44.48   male     45.05   male     46.98   male      The code:               # the GLM     > MF2 <- glm(Sex ~ Tarsus, data=mus, family=binomial)     > summary(MF2)          #attempt to plot sex (male) vs probability of tarsus being a specific length     > D1<-data.frame(Tarsus= mus$Tarsus[mus$Sex==1])     > p1<- predict(MF2 newdata=D1)     Error: unexpected symbol in "p1<- predict(MF2 newdata"     > p1<- predict(MF2, newdata=D1)     > lines(D1$mus, p1, lty=1)     > predict(MF2)


title: Why is a zero-intercept linear regression model predicts better than a model with an intercept?

text: Many textbooks and papers said that intercept should not be suppressed. Recently, I used a training dataset to build a linear regression model with or without an intercept. I was surprised to find that the model without an intercept predicts better than that with an intercept in terms of rmse in an independent validation dataset. Is the prediction accuracy one of the reasons that I should use zero-intercept models?


title: How can I analyse results from repeated measures of mixed design in R

text: I want to investigate performance in a learning test, where my test animals have been submitted to a number of trials within test day over several test days. I want to investigate if there is an effect of trial and of test day on the performance. Normally, I would use a random intercept model using `lme` in `R` with test animal as random factor, but I have not been able to find a method to analyse the effect of the two above-mentioned time factors in the same model. I would be grateful for anything pointing me in the right direction.


title: For what distribution is a trimmed mean the maximum likelihood estimator?

text: The sample mean is the maximum likelihood estimator of $\mu$ for a normal distribution $\text{Normal}(\mu,\sigma)$. The sample median is the maximum likelihood estimator of $m$ for a Laplace distribution $\text{Laplace}(m,s)$ (also called the double exponential distribution). **Does a distribution exist with a location parameter that the trimmed sample mean is the maximum likelihood estimator for?**


title: Are mathematically-coupled variables spurious predictors in linear regressions?

text: The correlation between two ratios with the same denominators are spurious. Similarly, the correlation between two mathematically-coupled variables could also be spurious. Are mathematically-coupled variables (e.g. age and income are used to derive two formulas, which are then used as Y and X in linear regression analyses) spurious predictors in linear regressions?


title: Gretl - calculate the estimated variance of the residuals

text: I want to learn the statistics package gretl. My first attempt to do so is to calculate a linear regression model of a set of data: $$y_i = \alpha + \beta x_i + u_i$$ First I want to create a crossplot of the data and then calculate the variance of the residuals $s^2_u$ and also those of $\alpha$ and $\beta$. My questions are: How to get the crossplot, I do not see any option in the menu which says crossplot? Is this possible in gretl? When I click Summary statistics I think I get the $s^2$ of the data, however is it also possible to get the the variance of the residuals $s^2_u$ and also those of $\alpha$ and $\beta$? Is it also possible to calculate in `gretl` a model without $\alpha$, such as: $$y_i = \beta x_i + u_i$$ I appreciate your answers! **UPDATE** My gretl window: ![enter image description here](http://i.stack.imgur.com/OSRkv.png)


title: Can someone explain this uncertainty estimation technique to me

text: I have been trying to understand a method of finding parameter uncertainties resulting from maximum likelihood estimation. Unfortunately the document I have is not in the public domain however I have found what I think is exactly the same method used in a paper from medicine. Having found it somewhere else has made me more confident that the method is correct but I still don't follow how/why it works. The paper is on pub-med at http://www.ncbi.nlm.nih.gov/pubmed/11316558 The section I don't get is 'calculation of the probability distribution of the parameters', on page 150. Is this a way of creating a 'credible interval' assuming a uniform prior? I have never dipped into credible intervals so I don't really follow. It's the 'normalising the parameter space to 1' that confuses/worries me. The section I don't understand proceeds as follows. 1) Calculate the likelihood L(x| a,b) for each point on the grid of parameter a vs. parameter b - so far so good. The next step is where I get lost, the calculation of P(a,b|x) from the above grid. 2) Sum together all the values on the grid and divide the value in each cell by this total - producing what I think is a 2D histogram with integral 1. The paper asserts that as P(a,b|x) is proportional to L(x|a,b) - which I agree with - then this normalised grid gives P(a,b|x). Is this simply a 'credible interval' with a uniform prior of 1. Perhaps what I don't understand is how the denominator in Bayes' theorem P(x) is given by the integral of a grid of L(a,b|x) for a range of a and b. Surely just by moving the grid one would change the values.


title: what is meant by stochastic drift?

text: could anybody define what is meant by stochastic drift? I _think_ I have a rough idea, say on a random walk X(n)= a + X(n-1) + Z(n) where Z(n) are iid zero mean and constant variance, then E[X(n)] = na. Is this stochastic drift? Many thanks


title: Seeking a new note-taking tool, for studying math and stats (with Latex support)

text: When studying math / stats or reading a academic paper often there are theorems, concepts, hypotheses and ideas are inter-connected to each other. These relationships can be a combination of hierarchical and flat network ones. From my own experience the best way to study is to taking notes. So, I am seeking a note-taking tool which I can 'map' the structure and contents while reading a paper/book. A similar tool I used for Mind Mapping is "Mindmeister", but it does not support Latex, and it's not very good at handling complex relationship as the requirement I described above. Does anyone have any suggestions? Thanks.


title: Independent and dependent variables in Repeated Measures ANOVA

text: I am working on a study where I have a sample of 86 students. I have tests for 5th grade, 6th grade, 7th grade, and 8th grade. I think I can run a Repeated Measures ANOVA since I have observations for 4 different occasions. However, I am unsure what my independent and dependent variables are in this set up.


title: Sign of correlation of logged variables

text: If we know that random variables $X$ and $Y$ are positively correlated, do we know that the random variables $Z = \log(X)$ and $W = \log(Y)$ are positively correlated?


title: Kolmogorov-Smirnov with discrete data: What is proper use of dgof::ks.test in R?

text: Beginner questions: I want to test whether two discrete data sets come from the same distribution. A Kolmogorov-Smirnov test was suggested to me. Conover ( _Practical Nonparametric Statistics_ , 3d) seems to say that the Kolmogorov-Smirnov Test can be used for this purpose, but its behavior "conservative" with discrete distributions, and I'm not sure what that means here. DavidR's comment on another question says "... You can still make a level α test based on the K-S statistic, but you'll have to find some other method to get the critical value, e.g. by simulation." The version of ks.test() in the dgof R package (article, cran) adds some capabilities not present in the default version of ks.test() in the stats package. Among other things, dgof::ks.test includes this parameter: > simulate.p.value: a logical indicating whether to compute p-values by Monte > Carlo simulation, for discrete goodness-of-fit tests only. Is the purpose of simulate.p.value=T to accomplish what DavidR suggests? Even if it is, I'm not sure whether I can really use dgof::ks.test for a two- sample test. It looks like it only provides a two-sample test for a continuous distribution: > If y is numeric, a two-sample test of the null hypothesis that x and y were > drawn from the same continuous distribution is performed. > > Alternatively, y can be a character string naming a continuous (cumulative) > distribution function (or such a function), or an ecdf function (or object > of class stepfun) giving a discrete distribution. In these cases, a one- > sample test is carried out of the null that the distribution function which > generated x is distribution y .... (Background details: Strictly speaking, my underlying distributions are continuous, but the data tend to lie very near to a handful of points. Each point is the result of a simulation, and is a mean of 10 or 20 real numbers between -1 and 1. By the end of the simulation, those numbers are nearly always very close to .9 or -.9. Thus the means cluster around a few values, and I am treating them as discrete. The simulation is complex, and I have no reason to think that the data follow a well-known distribution.) Advice? Thanks!


title: How to do PCA with a set-valued variable?

text: I wish to do a principal component analysis (PCA) with a few features (variables). While most of the features have a single numerical value, one feature is an _unordered_ vector of 5 numbers. For example:               #observation   value_feature1   value_feature2   value_feature_3     obsv1           1.5             0.65             (1.8,2.5,3.0,2.1,2.2)     obsv1           2.0             0.72             (2.1,2.2,1.8,1.7,2.5)      Is there any way I can do PCA with this type of data? Or do I need to combine the values for `feature_3` into one single value by some way; e.g., by taking the sum or product etc.? But in that case it might lose its physical meaning. In a second case `feature_3` can be a set of numbers but with variable set size, e.g. (1.2, 1.6) or (2.1, 2.3, 1.9) etc. How to do a PCA in that case?


title: density estimation when only mean and variance are known

text: I want to estimate the distribution function of the random variable. Due to the complicate nature of this function, I can be calculated only the first and second order moments(mean and variance) of this variable. Given these circumstances, how can I estimate the approximate distribution function? Thanks to all who help me to solve this problem, I am waiting for your helpful answers.


title: Cronbach's alpha negative result

text: What should I do if I am getting negative results for Cronbach's alpha? Is there anyway to adjust it?


title: How can posterior be persisted and reconstituted as future prior?

text: Suppose I model a data generating process as a hierarchal model and have made some training observation from the process. To learn about the process, with the observations I run the bayesian machinery to calculate posterior for the stochastic parameters of the model. Now I would like to be able to throw out the observations and just use the posterior I got in the learning process as prior for learning from the next set of observations. For cases where I don't have conjugate prior story, I am not sure how one can persist the state of posterior (while throwing away all your observations) so that it can be used as prior for the next observation? Any thoughts?


title: Negative inner fence value in a box and whiskers plot?

text: I calculated a negative lower inner fence value. I've double checked my work but I can't figure out what's wrong. I was under the impression that these values couldn't be outside of the min and max (duh). My values range from 0.06188 to 1.90957 and when I plot them in Graphpad it looks fine. I'm plotting with the 10% and 90% percentile whiskers (calculated as 0.515 and 1.359 respectively by Graphpad). I'm calculating the fraction of my datasets that are outside of the Tukey range (+/- 1.5*IQR), and the 10% and 90% range (possibly outliers). But when I calculate the Tukey fences myself using               lower inner = Q1-1.5(IQR)     upper inner = Q3+1.5(IQR)      I get a negative number for the lower value (-0.06758). I calculated the 10%/90% using Graphpad, so that's not a problem for now. Here's are the values, if anyone wants to recreate my calculations.               0.06187564,0.2931995,0.4433044,0.4908953,0.4992562,0.5123487,0.5188922,0.5326462,     0.5400167,0.5722309,0.6232626,0.6340977,0.6361935,0.6417624,0.6618853,0.6727413,     0.6990075,0.7250414,0.7575482,0.8221549,0.8573231,0.8732418,0.8749654,0.8751485,     0.88039,0.8831163,0.8842309,0.8927583,0.9143585,0.9371856,0.9644587,0.9662355,     0.9758021,0.9919777,1.004403,1.027557,1.037168,1.057354,1.075599,1.110754,1.120563,     1.135613,1.135632,1.135986,1.137239,1.148205,1.156857,1.175097,1.197407,1.21761,     1.219516,1.265934,1.304216,1.310815,1.336853,1.341317,1.357579,1.359387,1.462214,     1.552785,1.621189,1.760872,1.909566      It has to be something dumb I'm missing.


title: Multiple Imputation for Mixed Effects models

text: I stumbled accross this related question from 2010, and I wonder if there has been any progress on using multiple imputation for mixed effects models ? I prefer to use R, though Stata is also available to me.


title: Logistic regression is slow

text: I am relatively new to machine learning and I have a data classification problem where each sample has ~1500 features (continuous) and the category is binary. I want to apply logistic regression, for learning purposes. I have done a vectorized implementation in Python, which works very well and is fast on other datasets with a smaller number of features (order of magnitude 10s-100s features). On my dataset, however, I can't even finish the training (on ~3000 samples), as it simply takes ages. Could this be a problem of: 1) my vectorization, 2) the number of features in a sample (time for some dimensionality reduction?), 3) the fact that the problem is not suited to logistic regression, 4) using Python, or something else altogether?


title: How to interpret the T-depiction of multiple comparisons used in R's multcompView package?

text: R's multcompView package offers a seemingly useful set of tools to visualize the results of multiple comparisons, e.g. Tukey HSD analysis, including boxplots and 'letters' to identify significantly-different groups. It also provides a "t" depiction, in which groups are labeled with either a box, a sort of pyramid, or nothing. For instance:               require(multcompView)     y <- c(0,1,0, 1,2,1, 2,3,2, 8,9,10)     df <- data.frame(l=gl(4,3), y=y)     multcompBoxplot(y~l, data=df)      I can't make any sense of these "t"s, and can't find a clear explanation of them in the literature or documentation. Can anyone explain them to me?


title: Classification vs. regression for prediction of the sign of a continuous response variable

text: Say I want to predict whether or not a project will be profitable. In my sample data, the response variable is actually a continuous variable: the $ profit/loss of the project. Because my ultimate goal is just a binary classification (profitable project or unprofitable project), should I use a classification technique? Or should I use regression so that I'm not throwing away the additional information a continuous response variable provides?


title: Separate points in 3 or more dimentions into 2 groups with labels A and B

text: How can you separate a group of points in 3 dimensions into two groups with labels A and B? My approach:   If we draw a plane through the cloud of points, the group A can be the points above the plane, and the group B can be the points below the plane. I am not sure if this a legitimate approach. Is anyone aware of any other methods, preferably non optimization, to do so?


title: How could I visualise the importance of different inputs to the forecast for a black-box non-linear model?

text: I am building an interactive forecast tool (in python) as an aid to forecasting that is done in my organisation. To date the forecast process has been largely human driven, with forecasters assimilating the data in their natural neural networks and using their learned gut feel to make predictions. From a long term forecast verification and predictive modelling study I've done I've found what you might expect; different forecasters exhibit different biases, the effects of some predictors seem to be overstated and other important ones seem to be ignored and in general the forecast performance is mediocre compared with relatively simple empirical models. The forecasts will continue to be manual, but I am trying to build a useful tool to provide the forecasters with a better quantification of the relative effects of predictors. There are also important effects such as seasonal influences that are often overlooked that I would like the tool to highlight to the user. I am expecting a degree of backlash and scepticism about the modelling process from some of the more 'experienced' forecasters (many of whom have little formal knowledge of statistics), so the communication is at least as important and the model performance itself in terms of achieving a measurable improvement in forecast accuracy. The models I'm developing have a strong auto-regressive component that is at times modified significantly by events which show up as measured values in some predictors that are, during non-event times, close to zero. This accords with the mental model that forecasters use. The key part is being able to demonstrate which of the 'event' measurements are most influential in driving the prediction away from the auto-regressive value for any given forecast. I imaging the process in this way; the forecaster divines their best guess value, the model suggests a different one and the forecaster asks why. The model replies something like "see here, this value of this predictor increases the forecast value in Summer. If it was Winter, it would move the other way. I know there are these other measurements, but they have much less effect than this one". Now, imagine the model was a simple linear regression. One could imagine displaying the relative 'effect' of event based predictors by multiplying the value by the model co-efficient and displaying as a simple bar chart. All the bars from the different predictors add up to the total deviation from the AR value, and this succinctly and clearly shows the ones that are, in this instance, having a strong influence. The problem is that the process being forecast displays a high degree of non- linearity in the predictors, or at least, I have had much more success with black-box non-linear machine learning algorithms (random forest and GBM) than with GLMs for this data-set. Ideally I would like to be able to seamlessly change the model working 'under the hood' without the user experience changing, so I need some generic way of demonstrating in a simple fashion the importance of the different measurements without using some algorithm specific approach. My current approach will be to quasi-linearise the effects by setting all values to zero except for one predictor, record the predicted deviation and then repeat for all predictors, displaying the results in the bar chart mentioned above. In the presence of strong non-linearity, this may not work so well. Are there any known approaches for achieving in a clear way what I am trying to do here?


title: Goodness of fit test: question about Anderson–Darling test and Cramér–von Mises criterion

text: I'm reading web pages for goodness of fit tests, when I came to the Anderson- Darling test and the Cramér-von Mises criterion. So far I got the point; it seems the Anderson–Darling test and the Cramér-von Mises criterion are similar, just based on a different weighting function $w$. Also there's a variant of the Cramér-von Mises criterion named the Watson test. Basically I have two questions here   1. There are not many Google results about these two methods; are they still state-of-the-art? or replaced by some better approaches already? It's a bit of a surprise, as according to this paper on power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests, AD is performing quite well; always better than Lilliefors and KS, and very close to the SW test, which is specifically designed for the normal distribution.   2. What is the confidence interval for such tests?  For the AD, CM and Watson tests, I saw the test statistics variable defined on the wiki pages, but didn't find the confidence interval. Things are just more straightforward for the KS test: on the wiki page, the confidence interval is defined by $K_\alpha$, which is defined from the cumulative distribution function of $K$.


title: Mining cycles from graphs

text: Suppose $T$ is a set of transactions containing some paths taken from a set of undirected graphs $G$ whose nodes are all taken from the same finite alphabet $A$. I would now like to extract all frequent cycles (strongly connected components) from $T$ by using the Apriori algorithm. The idea is that since the Apriori is anti-monotone in the support parameter (for the item $(abc)$ to be frequent, all subsets $(ab)$, $(ac)$, $(bc)$ etc must also be frequent), it is sufficient to just extract all rules of length $2$, ie. $a\to b$, $b\to c$ etc. Returning to the graph representation, we do not need to extract paths like $a\to b \to c \to d$, because they can all be described as a set of paths of length $2$. Am I correct in saying that by combining all such rules/paths of length $2$ from $T$ into a new directed graph $G'$, the strongly connected components can be extracted?


title: Using continuity correction for normal approximation or not?

text: Below is a question on a recent actuarial exam, Exam 3L of the CAS. I didn't know whether or not to use the continuity correction when using the normal approximation to do hypothesis testing involving a bernoulli trial. The answer is vastly different (based on the answer choices given in the exam) depending on whether or not you use it. The problem does not specify to use it or to not use it, though some previous problems have specified one or the other, whereas others have not specified either. > You are given the following: > >   * Accidents happen during a work day at a probability of p when a machine > is operated. >   * The null hypothesis $H_0$ is that the probability of an accident is > 0.05; the alternative hypothesis $H_1$ is that the probability is less than > 0.05. >   * If less than 20 accidents are observed in 365 work days, then reject the > null hypothesis. > > > Using the normal approximation, calculate the probability of Type II error > using the value 0.03 as the true probability of an accident occurring. Any help would be greatly appreciated. Just to be clear, I am hoping to better understand when it should be used and when it should not be used, in general, in addition to learning the best way to do this problem.


title: User-defined correlation matrix in R package nlme with negative values

text: I have a nonlinear model with residuals that are negatively autocorrelated at short distances. I can find no spatial correlation structures in `nlme` that can easily handle negative autocorrelation as most have bounds on parameter values so that the correlation is between 0 and 1. First, is there something I am missing? I tried to roll my own by calculating a correlation matrix with some negative values off the diagonal and then setting the correlation structure as follows within the function `gnls`:               correlation=corSymm(corr9x[lower.tri(corr9x)])       where `corr9x` is the $n \times n$ matrix of correlations that I set up based upon distances between points in the data set. Some of these correlations are positive and some are negative. They are based on a Moran's I correlogram that I calculated from the residuals returned from a `gnls` model fit with NO spatial correlation employed. I get the following error:               Initial values for corSymm do not define a positive-definite correlation structure       I am unsure if the matrix is rejected out of hand because it contains negative values or if there is something I can do to coerce it. I have checked the lower triangle matrix returned and it matches what I intended. Any input is appreciated.


title: Standard deviation of error in simple linear regression

text: How to derive the following formula: where $\sigma(y)$ is the standard deviation of $y$ (dependent variable), and $\rho(x,y)^2$ is correlation between $x$ and $y$ squared, $\sigma(\epsilon)$ the standard deviation of the error. $$ \sigma(\epsilon) = \sigma(y) \sqrt{1-\rho(x,y)^2} $$ if $y = a +bx + \epsilon$.


title: Simulating Monte Carlo with different standard deviations and interval confidence

text: I have a question regarding Monte Carlo simulation (direct simulation), applied to propagation of uncertainties. From what I understand Monte Carlo accepts random numbers of each input variable of the model. These random numbers are generated with the mean, standard deviation and type of PDF (normal, triangular, uniform, etc.) But the question is: I have 2 input variables in my model:               x: mean=30.5 std=0.001 <- calculated with 95% confidence.     y: mean=2 std=0.1 <- calculated with 99% confidence          with z = x*y      What confidence interval will `z` be generated? How do I know which uncertainty (standard deviation) of `z` is correspondent to a 95% confidence level?


title: Kernel bandwidth in Kernel density estimation

text: I am doing some Kernel density estimation, with a weighted points set (ie., each sample has a weight which is not necessary one), in N dimensions. Also, these samples are just in a metric space (ie., we can define a distance between them) but nothing else. For example, we cannot determine the mean of the sample points, nor the standard deviation, nor scale one variable compared to another. The Kernel is just affected by this distance, and the weight of each sample: f(x) = 1./(sum_weights) * sum(weight_i/h * Kernel(distance(x,x_i)/h)) In this context, I am trying to find a robust estimation for the kernel bandwidth 'h', possibly spatially varying, and preferably which gives an exact reconstruction on the training dataset x_i. If necessary, we could assume that the function is relatively smooth. I tried using the distance to the first or second nearest neighbor but it gives quite bad results. I tried with leave-one-out optimization, but I have difficulties finding a good measure to optimize for in this context in N-d, so it finds very bad estimates, especially for the training samples themselves. I cannot use the greedy estimate based on the normal assumption since I cannot compute the standard deviation. I found references using covariance matrices to get anisotropic kernels, but again, it wouldn't hold in this space... Someone has an idea or a reference ? Thank you very much in advance!


title: Density estimation methods?

text: What non-/semiparametric methods to estimate a probability density from a data sample are you using ? (Please do not include more than one method per answer)


title: Meaning of "a statistic $U$ is ancillary to another statistic $T$"?

text: From Wikipedia > Given a statistic $T$ that is not sufficient, an ancillary complement is a > statistic $U$ that is **ancillary to $T$** and such that $(T, U)$ is > sufficient. Intuitively, an ancillary complement "adds the missing > information" (without duplicating any). I understand that an ancillary statistic wrt a family of possible distributions of the sample, is defined as a statistic whose distribution doesn't depend on any sample distribution in the family of sample distributions. I was wondering what "a statistic $U$ is ancillary to another statistic $T$" mean, in the definiiton of an ancillary complement? Thanks and regards!


title: Can you use normal correlation for vectors with only 2 (or 3), ordered, levels?

text: When one wants to compute the correlation of two vectors of a continues variables, one uses pearson (or spearman) correlation. **But what should (can) one use for the case of two vectors with 2 (or 3) ordered levels only? Is spearman enough, or does it require another method?** I remember coming across someone who once claimed to me that OR (odds ratio) is more fitting for such situations (for 2 by 2 tables, where order has no meaning), is this true? Here is an example R code, for allowing of answers relating to the same example:               set.seed(10)          x2 <- sample(c(-1,1),50, T)     x3 <- sample(c(-1:1),50, T)     y3 <- sample(c(-1:1),50, T)     y2 <- sample(c(-1,1),50, T)          cor(x3,y3, method = c("spearman"))     cor(x2,y2, method = c("spearman"))     cor(x3,y2, method = c("spearman"))      p.s: for the 2 by 2 case, I followed from the comments that categorical "measures of association" is the term to look for. However, Part of the time I am comparing 2 on 3 tables, where on the factor with 3 levels there is order - so I would like to take use of that information.


title: Question on autoregressive equation

text: I was given a reduced form VAR model, where the dependent variable is inflation, and independent variables include inflation lagged by four periods (L.Inflation) and other exogenous variables. The lagged part is the AR component of the model. Shouldn't the AR include L1.Inflation, L2.Inflation, L3.Inflation as well? Is this still an AR equation without these three terms?


title: Books covering how to conduct a meta-anlysis

text: I would like to know if anyone could recommend a book that deals more with the practical issues around conducting a meta-analysis? Thanking you in advance Andrew Vitiello


title: Countering false positives in GWAS due to linkage disequilibrium

text: We know that due to LD, we can get significant p-values for markers near a causal marker (or a marker closes to the causal region) in GWAS studies. I've seen attempts looking at LD to do multiple testing correction, by counting the effective number of tests, but as far as I know, they are simply trying to fix the correct global p-value cutoff for a certain error rate. My question is regarding the actual problem of selecting causal markers. What attempts have their been to remove these false positives that occur due to them being in LD with the causal markers? It seems regression and penalized regression will remove some correlated markers, but are their any techniques other than penalized regression to remove these false positives? Regression is obviously also intractable when we have a large number of markers.


title: The exact distribution of Wilcoxon rank-sum statistic U

text: The distribition of the rank-sum statistic U is assumed to be normal for large number of samples being considered. What is the exact distribution? I want to compare and sometimes fuse results from various tests wherein some tests might not have large number of samples. I want to have a the exact distributions in cases where, say $n_1 n_2 < 30$. Is there a closed form that can be used or calculated? **Update** : So apparently, people cite Streitberg, B. and J. Rohmel, Exact distributions for permutation and rank tests: An introduction to some recently published algorithms, Statist. Software Newsletter 1 (1986) 10-17. for the exact distribution, but I have neither been able to find it nor the result yet.


How can I generate events using the Poisson distribution in R? The events could be the occurrence of floods in the next 1000 years at a given rate of occurrence per year.


title: SVM with quadratic loss

text: I've seen some statement where I got the impression that SVM with a quadratic loss is no more than having a kernel matrix where a multiple of the unit matrix is subtracted from the kernel. It was shown that in the dual problem the L2 loss constant can be merged into the kernel. But it wasn't clear to me. So the question: Are L1 and L2 loss for SVM substantially different? (Is L2 loss for SVM convex?) Or does L2 loss correspond to L1 loss with a shifted kernel? Or anything else?


title: Necessary amount of countries in the model for 1 macro level predictor?

text: Within my logistic regression models I use cross-national data in order to say something about individual soft drug tolerance. In total I have 29 countries with 37.000 individuals. My supervisor wants me to find some theoretical fundament to the rule of thumb I have used (and she proposed), namely the fact that for every 10/15 countries 1 possible macro-level predictor can be included in the regression model. Can anybody provide me with this literature (and a reaction to this presumed rule of thumb). Thanks in advance!


title: References about the theory of linear regression or regression in general

text: I am looking for references about the theory of linear regression or regression in general. More specifically, I am interested in knowing under what circumstances an estimated regressor is going to give accurate solutions with certain probabiltiy. For now, I am looking for the basic classical results, nothing too fancy or modern.


title: Search for interactions using carets rfe function

text: it is straight forward to search for purely additive models using the rfe function in caret. Is it possible to include all interactions as part of the search? In the train method, we can simply say train(Y~.^2, data) to include interactions, but I have not found a way to do this in rfe. I could probably just make a data frame where all interactions are included and use that in rfe, but I'm hoping that there is a more elegant solution. Thanks.


title: Which bootstrapped regression model should I choose?

text: I have a binary logistic regression model with a DV (disease: yes/no) and 5 predictors (demographics [age, gender, tobacco smoking (yes/no)], a medical index (ordinal) and one random treatment [yes/no]). I have also modeled all the two-sided interaction terms. The main variables are centered and there is no sign of multicollinearity (all VIFs < 2.5). I have some questions:   1. Is bootstrapping advantageous over my single model? if so,   2. which bootstrapped model should I choose? I just wanted to see if bootstrapping algorithms follow random methods for creating new samples, or if they have rigid algorithms. Therefore, I resampled for 1000 times in each attempt (so I have several bootstrapped models, each with 1000 trials). However, each time the coefficients of the bootstrapped model differ (although the number of trials are constantly 1000). So I wonder which one should I choose for my report? Some changes are tiny and don't affect my coefficients' significance, but some make some of my coefficients non-significant (only those with P values close to 0.05 in the original model that change to 0.06 for example).   3. Should I choose a higher number like 10,000? How can I determine this limit?   4. Again should I bootstrap in the first place? If its results vary each time, can I rely on its results?   5. Do you have any other ideas in mind that can help me with my case? Many many thanks.


title: Determining Optimal Number of Cluster in Hierarchical Clustering in Consideration of Variance of Data

text: I'm applying a Hierarchical Agglomerative Clustering (HAC) for grouping my data and I need to determine the number of the cluster automatically. To determine the optimal number of cluster, I obtain the the best cluster combination which maximizes the similarity of each member in one cluster and minimizes the similarity between clusters. So far my strategy is working well. But I have one case where my HAC algorithm is totally failed, as you can see in this image there are two groups (in upper right) which the variance of each members is really small (very small compared to others...). I expect my clustering would end like this, but unfortunately, this is what i get. The data with larger variance, ends up with small partitions instead group them into three bigger clusters because of these two groups, which each data is really similar (almost the same) with each other, it messed up with optimization strategy to determine the best number of cluster of my HAC algorithm. I want to ask if there's any optimization strategy especially for HAC which can deal with this kind of case. So far I'm still looking for any strategy which considers variance of the data but I'm still not sure, Currently the best strategy that popped out in my head is :   1. Automatically detect the sub-clusters (or sub-tree of dendogram of HAC) which are considered having small variance (or the similarity is really close)    2. Discard them    3. Recalculate HAC.  The problem is how do I define "small value of variance"? If you know any better strategy i would be happy to know :) Thanks


title: Is this a "statistical matching" problem?

text: For a finite population I have two data sources A and B with variables (X, Y) and (X, Z), respectively. The first data source is actually a full sample, the second is a simple random sample. The items in B cannot be identified precisely using X alone. I want to combine these data sets using the common variable X, assuming conditional independence of Y and Z given X in the simple case and perhaps later using auxiliary information on the interdependence of X, Y and Z. The resulting data should be an extension of the full sample A.               +---+---+     |   B   |     +---+---+---+         |   A   |     +---+---+---+       +---+---+---+     | Z | X | Y |       | X | Y | Z |     +---+---+---+       +---+---+---+     |   | x | A |       | x | A | B'|     | B | x | A |       | x | A | B'|     |   | x | A |       | x | A | B'|     |   | x | A | ----> | x | A | B'|     |   | x | A |       | x | A | B'|     | B | x | A |       | x | A | B'|     |   | x | A |       | x | A | B'|     +---+---+---+       +---+---+---+      From the literature that treats data fusion and statistical matching (D'Orazio et al. 2006, Kiesl and Rässler 2006) I understand that my problem does not fall in this category, since A and B are not disjoint. Kiesl and Rässler call this "database enrichment" (case 2 in Figure 1, page 4), but a search for this term didn't find anything useful. Is there a more specific name for the task I am facing? Could you point me at literature?


title: Judging Matthews correlation coefficient

text: I often read the Matthews correlation coefficient is one of the most generally used validity estimates of confusion matrices in machine learning. However, I couldn't find a reference which states what a good value of this coefficient should be. What is the status on that issue?


title: Correlation of variable corrected for age in SPSS

text: We have measured reading speed in patients with diabetic eye disease. What we have found is that reading speed is faster in younger patients and in those who has better vision as you would expected. However, we have developed an novel measurement which is much faster to do than measuring reading speed. Both reading speed and the novel measurement correlate well with each other but we would like to see whether this correlate is independent of age and vision level. How can I do that in SPSS. Many thanks.


title: Tutorials / examples for multiclass boosting

text: I want to learn the multiclass boosting technique. I have a basic understanding of binary boosting and also have seen some working examples on this. I have also read about the basics of multiclass boosting algorithms. But I have not been able to search for any examples in multiclass boosting. Can anyone please guide me How to learn on this topic?