title: How to calculate Standard Error of Odds Ratios?

text: I have two datasets from genome-wide association studies. The only information available is the odds ratio and the p-value for the first data set. For the second data set I have the Odds Ratio, p-value and allele frequencies (AFD= disease, AFC= controls) (e.g: 0.321). I'm trying to do a meta-analysis of these data but I don't have the effect size parameter to perform this. Is there a possibility to calculate the SE and OR confidence intervals for each of these data only using the info that is provided??   Thank you in advance example: Data available:                   Study     SNP ID      P        OR    Allele   AFD    AFC         1         rs12345    0.023    0.85         2         rs12345    0.014    0.91     C      0.32   0.25      With these data can I calculate the SE and CI95% OR ? Thanks


title: MATLAB interpretation of xcorr2 cross correlation function

text: I have two vectors of matching lengths. They are readings from two different sensors (one is from a smartphone and the other is from a wiimote) of the same hand movement. I am trying to find the time offset between them to synchronise the readings for further processing. The readings I get are of the format (Time(ms) Value) for accelerations in the X,Y and Z direction. For the synchronization, I plotted the cross-correlation function `xcorr2()` between the two sets. I am getting the same graph (a weird triangle peak and a straight line at the bottom) for Accelerations along the x, y and z directions (which I guess is good) but I don't know how to interpret it. What do the axes in the graph represent? Can anyone explain to me what `xcorr2()` means in a qualitative sense. From the correlation function, how do I determine the offset (i.e. how many seconds is sensor1 behind sensor2)? Thanks! Imelza


title: JAGS: posterior predictive distribution

text: I am fitting a simple linear regression model with R & JAGS:               model {         for (i in 1:length(ri)){             ri[i] ~ dnorm(mu[i],tau)             mu[i] <- alpha + b.vr*vr[i] + b.ir*ir[i]         }              #posterior predictive distribution         for (j in 1:length(pvr)){             pmu[j] <- alpha + b.vr*pvr[j] + b.ir*pir[j]             pri[j] ~ dnorm(pmu[j],tau)          }                #priors for regression         alpha ~ dnorm(0,0.01)         b.vr ~ dunif(-10,0)         b.ir ~ dunif(0,10)              #prior for precision         tau ~ dgamma(0.001,0.001)         sigma <- 1/sqrt(tau)     }      I am trying to calculate the posterior predictive distribution with new data (`pvr` and `pir` in the model). But somehow the reported results (`pri`) do not make sense (the means of `pri` are smaller than expected). Could someone explain me, is something wrong with the model specification?


title: Algebraic Geometry for statistics

text: I have heard about uses of Algebraic Geometry in statistics and Machine Learning. I wanted to try to learn a bit about this topics. I don't know nearly anything about algebraic geometry, but I have background in math, and I know about basic group theory, rings fields and some commutative algebra. My questions are:   1. What are the Algebriac Geometric concepts I should learn that are connected to applications in stat/ML (I suppose only a portion of what is usually taught in Algebraic Geometric courses and books is useful).   2. Can you recommend some books / introductory papers for someone with my background? I don't mean standard textbooks for AG but something that focuses on concepts used in applications.


title: Matlab - nu-support vector regression with LibSVM

text: I've 2 problems when applying $\nu$-SVR with LibSVM in Matlab.   1. How to extract support vector from training model in $\nu$-SVR?   2. When I visualize `slab tube` for a training and/or a testing data. I take $\varepsilon$ value in model and plot $f(x) - \varepsilon$ and $f(x) + \varepsilon$ for upper and lower parts. But, according to `Scholkopf (2000)`:  `Let us briefly discuss how ν-SVR relates to ε-SVR (see section 1). Both algorithms use the ε-insensitive loss function, but ν-SVR automatically computes ε.` I've check the $\varepsilon$ value before and after training, but it's not changed. So, how is LibSVM applying this `ν-SVR automatically computes ε`?


title: Interpretation of Spearman correlation for small sample

text: I am currently being confused by different opinions regarding the spearman correlation interpretation. Some says, $\rho<0.2$ can be ignored as small/weak relationship, some says $\rho<0.1$, for sample size. I am explaining it here with examples. Suppose, I have 29 samples with $x$ and $y$ features, $\rho_{x,y}=0.533$. I did two transformation on $y$ to make this relationship $<0.1$ and compare both transformation. This is to remove the influence (relationship) of $y$ with $x$. What I have obtained is 1) Transformation 1: $\rho_{y_{t1},x}=0.54$ 2) Transformation 2: $\rho_{y_{t2},x}=-0.13$ For this small sample size can I safely conclude that transformation 2 successfully reduced the influences of $x$ on $y_{t2}$ (to make it almost independent of $x$)? Thanks, Wahid


title: Subsetting R time series vectors

text: I have a time series and I want to subset it while keeping it as a time series, preserving the start, end, and frequency. For example, let's say I have a time series:               > qs <- ts(101:110, start=c(2009, 2), frequency=4)     > qs          Qtr1 Qtr2 Qtr3 Qtr4     2009       101  102  103     2010  104  105  106  107     2011  108  109  110           Now I will subset it:               > qs[time(qs) >= 2010 & time(qs) < 2011]     [1] 104 105 106 107      Notice that I got the correct results, but I lost the "wrappings" from the time series (namely start, end, frequency.) I'm looking for a function for this. Isn't subsetting a time series is a common scenario? Since I haven't found one yet, here is a function I wrote:               subset.ts <- function(data, start, end) {       ks <- which(time(data) >= start & time(data) < end)       vec <- data[ks]       ts(vec, start=start(data) + c(0, ks[1] - 1), frequency=frequency(data))     }      I'd like to hear about improvements or cleaner ways to do this. In particular, I don't like the way I'm hard-coding start and end. I'd rather let the user specify an arbitrary boolean condition.


title: Visualizing the intersections of many sets

text: Is there a visualization model that is good for showing the intersection overlap of many sets? I am thinking something like Venn diagrams but that somehow might lend itself better to a larger number of sets such as 10 or more. Wikipedia does show some higher set Venn diagrams but even the 4 set diagrams are a lot to take in. My guess as to the final result of the data would be that many of the sets won't overlap so it is possible that Venn diagrams would be fine -- but I would like to find a computer tool that will be able to generate that. Its looks to me like Google charts doesn't allow that many sets.


title: Which Test do I use? ANCOVA, repeated measures, multi-way ANOVA?

text: We are analysing the effects of different harvesting intensities on soil nutrients. Data was collected over three different years. The first year was before treatments were applied, and the second and third years were 1 year and 6 years post treatment. There are 4 stand types (each replicated 3 times), 6 treatments (retention harvesting of all different intensities) and 6 plots within each stand type*treatment. I am having some trouble wrapping my head around this question and am wondering which test I should be using?


title: How to prove $(X_{n})_{n\in \mathbb N}$ and $(Y_{n})_{n\in \mathbb N}$ are supermartingale and $(Y_{n})_{n\in \mathbb N}$ is convergence to -7

text: Let $p \in [0 , \frac{1}{2}] $ and $\eta_{i}$ be i.i.d random variables and $P(\eta_{i}=1)=p$ and $P(\eta_{i}=-1)=1-p$ and $\mathcal F_{n}=\sigma(\eta_{1},\cdots,\eta_{n})$ and $X_{n}=\sum_{i=1}^{n}\eta_{i}$ and $Y_{n}=X_{T(-7) \wedge n} $(means martingale/submartingale/supermartingale when stopped by a stopping time) . show that   1. $(X_{n})_{n\in \mathbb N}$ and $(Y_{n})_{n\in \mathbb N}$ are supermartingale or matigale or submartingale.


title: Contingency tables: what tests to do and when?

text: I'd like to see an extension of this discussion of the age-old chi-sq vs. Fisher's exact test debate, broadening the scope a bit. There are many many tests for interactions in a contingency table, enough to make my head spin. I'm hoping to get an explanation of what test I should use and when, and of course an explanation as to why one test should be preferred over another. My current problem is the classic $n \times m$ case, but answers regarding higher dimensionality are welcome, as are tips for implementing the various solutions in R, at least, in cases where it is non-obvious how to proceed. Below I've listed all the tests I'm aware of; I hope by exposing my errors they can be corrected.   * $\chi^2$. The old standby. There are three major options here:     * The correction built into R for 2x2 tables: "one half is subtracted from all $|O-E|$ differences." Should I always be doing this?     * "$N-1$" $\chi^2$ Test, not sure how to do this in R.     * Monte Carlo simulation. Is this always best? Why does R not give me df when I do this?   * Fisher's exact test.     * Traditionally advised when any cell is expected to be <4, but apparently some dispute this advice.     * Is the (usually false) assumption that the marginals are fixed really the biggest problem with this test?   * Barnard's exact test     * Another exact test, except I've never heard of it.   * Poisson regression     * One thing that always confuses me about glms is exactly how to do this significance tests so help on that would be appreciated. Is it best to do nested model comparison? What about a Wald test for a particular predictor?     * Should I really just always be doing Poisson regression? What's the practical difference between this and a $\chi^2$ test?


title: How is LSD an incomplete three-way layout?

text: How is the Latin square design(LSD) an incomplete three-way layout? In this context, let me ask also what is a complete layout?


title: Group-level dependent variable and within-group interdependent events (Mplus?)

text: Here is a minimum viable example that resembles my problem:   * Subjects are human individuals   * Subjects are characterized by a set of categorical and continuous variables   * Subjects engage in a varying number of consecutive events   * I assume that the number and kind of events they engage in depends to some extent on the subjects’ characteristics   * Each event is characterized by a set of categorical and continuous variables   * The dependent variable is "success", which is binary coded (1 = success / 0 = failure)    * The series of events for each individual ends either with success or failure.   * I assume that the characteristics of an event dependent to some extent on the characteristics of previous events and the individual characteristics.   * I assume that the event series outcome (1 / 0) is influenced by the subjects‘ characteristics, the single events’ characteristics and characteristics of the event series.  My goal: I am trying to develop a model that predicts the probability of success for a new individual at any given event in his /her "event series". The dependent variable is measured at group level (individual subject = group (or series) of events). Hence I am aware that I am facing a multi-level problem with a group-level dependent variable that needs to incorporate some form of Bayesian updating at the within-group level.   * The group-level variables act as predictors for group-level outcome (success) and the specific series of events (what are the events' characteristics, how are the events ordered, how many events are there)   * The within-group variables (the event characteristics) also act as predictors on subsequent event characteristics and the group-level outcome (success)


title: What kind of plot am I looking at?

text: I stumbled on to these following two slides (slides 21 & 22 on a machine learning tutorial found here): ![Scatterplot..](http://i.stack.imgur.com/0KPox.png) ![enter image description here](http://i.stack.imgur.com/8uoMa.png) The first is obviously an $x,y$ scatterplot of height and weight. **But what is the second?** Looks like two component parts representing the original data, but I can really tease out what is going on. Any insight (or reproducible R code, or even Pseudo-code) on what the heat map is representing, and what the white lines are interpreted as would be great.


title: Is $Max${$F(X),G(Y)$} necessarily a distribution function?

text: If $F(X)$ and $G(Y)$ denote two distribution functions, then is $Max${$F(X),G(Y)$} necessarily a distribution function as well?


title: Mean-reversion for random process

text: I moedelled following ARIMA Process which I also smoothed. > abc1 <- (arima.sim(n = 1400, list(ar = c(0.3, 0.2, 0.15, 0.1, 0.1, 0.1)), sd > = sqrt(0.5)))   >  plot(abc1)   >  abc2 <- SMA(abc1, n=360)   >  plot(abc2) My aim is it to add a component for mean-reversion to describe a business cycle model. Ideas? thanks, F!


title: How can I generate predictions from the randomSurvivalForest package in R?

text: I'm trying to use the `randomSurvivalForest` package in `R` to predict the next event in a series of events (using `res <\- predict(fit,v)`), but all I get is a total % survivors and then the cumulative hazard for each individual in the `res$ensemble`. Can I use this to predict/generate values matching what's in in `res$time`?


title: Confidence measure for classification result from a MAP estimator

text: I'm using a maximum a posteriori probability (MAP) estimator in a classification problem. After estimating all the a posteriori probability, the standard way is to simply take the class associated with the maximum probability. I would like to measure the confidence of this decision. I.e. if the maximum a posteriori probability is not significantly larger than the rest, I would like to skip the case (which would reduce recall but increase precision). Is there a standard way to measure such concept?


title: The "sum" of prediction intervals

text: Suppose I have a regression model which yields a couple of predicted values with their respective prediction intervals and the random quantity that I am interested in is the sum of (some subset) of those predicted values. How do I go about building a new prediction interval for this sum? Is it as easy as summing up the limits of the respective prediction intervals?


title: How to compare two or more correlation matrices?

text: I have $P$ correlation matrices $(n \times n)$ computed with $P$ sets of $(m \times n)$ data (observed) using the MATLAB function `corrcoef`.   * How do I compare and analyze these $P$ correlation matrices with respect to each other?    * What are the tests, methods and/or checkpoints?


title: What is the difference between functional data analysis and high dimensional data analysis

text: There are a lot of references in the statistic literature to " **functional data** " (i.e. data that are curves), and in parallel, to " **high dimensional data** " (i.e. when data are high dimensional vectors). My question is about the difference between the two type of data. When talking about applied statistic methodologies that apply in case 1 can be understood as a rephrasing of methodologies from case 2 through a projection into a finite dimensional subspace of a space of functions, it can be polynomes, splines, wavelet, Fourier, .... and will translate the functional problem into a finite dimensional vectorial problem (since in applied mathematic everything comes to be finite at some point). **My question is:** can we say that any statistical procedure that applies to functional data can also be applied (almost directly) to high dimension data and that any procedure dedicated to high dimensional data can be (almost directly) applied to functional data ? If the answer is no, can you illustrate ? EDIT/UPDATE with the help of Simon Byrne's answer:   * sparsity (S-sparse assumption, $l^p$ ball and weak $l^p$ ball for $p<1$) is used as a structural assumption in high dimensional statistical analysis.    * "smoothness" is used as a structural assumption in functional data analysis.  On the other hand, inverse Fourier transform and inverse wavelet transform are transforming sparcity into smoothness, and smoothness is transformed into sparcity by wavelet and fourier transform. This make the critical difference mentionned by Simon not so critical?


title: Sample size for LASSO

text: What considerations should be made regarding sample size for using LASSO or elasticnet? We are going to gather expression data from 1700 genes and our response variable is multinomial (3 categories).


title: Interpreting the Johansen cointegration test output

text: I am running the Johansen cointegration test using 2 non-stationary time- series, as suggested by the literature. The output I got is the following:                   r0        r1           t1  0         1      Particularly, i'm using the Matlab function jcitest() to run this test. In the link above, that explains in details how the function works, it says that the [h] output value corresponds to a Boolean decisions for the tests; particularly, values of h equal to 1 (true) indicate rejection of the null of cointegration rank r in favor of the alternative, while, values of h equal to 0 (false) indicate a failure to reject the null. My question now is: how can I establish if the time-series are cointegrated or not, since that with rank 0 they result to be cointegrated, while with rank 1 they are not? What is the rule to establish if the time-series are cointegrated or not? Thanks for help in advance.


title: Estimating successes while obtaining Bernoulli samples

text: I have a process which, after fixing the values of some parameters, generates samples from a Bernoulli distribution with unknown $p$. The value of $p$ is typically small, and what I want to do is to discover suitable values for my parameters so that $p$ is at least $0.1$. An added problem is that generating each sample (i.e. running my “experiment” once) takes some considerable amount of time. One thing I can do is, say, fix the parameters, generate 100 samples, count the number $k$ of successes and, if $k/100 < 0.1$, try again with different parameters until I find ones that yield $k/100 > 0.1$. However, as generating each sample takes some time, intuitively I would like to stop generating samples for a fixed set of parameters if it “doesn't look promising”. For example, say, if i've already seen 30 samples and not a single success then, according to this question, with 95% confidence the value of $p < 0.1$; so it would be reasonable to stop here and not generate the remaining 70 cases. I would like to generalise this idea. I guess the question I really want to ask is the following: > Given that I've already seen $k$ successes on $n$ samples of a Bernoulli > distribution with unknown parameter $p$, what is the probability that, if I > keep sampling from the same distribution, after $N$ observations (say $N = > 100$), I'll see at least $K$ successes (say $K = 10$). **Update** : In case anybody is interested, these are some example values I computed with a small script in python using the answers provided below. Note that both methods are trying to compute two _slightly_ different probabilities, you might want to read the exact details below.                     Rasmus' method                | Zen's method                            k=0  k=1  k=2  k=3  k=4  k=5  | k=0  k=1  k=2  k=3  k=4  k=5      n= 0  0.90                          | 0.90                              n= 1  0.81 0.99                     | 0.81 0.99                         n= 2  0.72 0.97 1.00                | 0.73 0.97 1.00                    n= 3  0.66 0.95 1.00 1.00           | 0.66 0.95 1.00 1.00               n= 4  0.59 0.92 0.99 1.00 1.00      | 0.60 0.92 0.99 1.00 1.00          n= 5  0.53 0.89 0.98 1.00 1.00 1.00 | 0.54 0.88 0.98 1.00 1.00 1.00     n= 6  0.48 0.85 0.97 1.00 1.00 1.00 | 0.49 0.85 0.97 1.00 1.00 1.00     n= 7  0.43 0.81 0.96 0.99 1.00 1.00 | 0.45 0.81 0.96 0.99 1.00 1.00     n= 8  0.39 0.78 0.95 0.99 1.00 1.00 | 0.41 0.78 0.94 0.99 1.00 1.00     n= 9  0.35 0.74 0.93 0.99 1.00 1.00 | 0.37 0.74 0.92 0.98 1.00 1.00     n=10  0.32 0.70 0.91 0.98 1.00 1.00 | 0.34 0.70 0.90 0.98 1.00 1.00     n=11  0.28 0.66 0.89 0.97 1.00 1.00 | 0.31 0.67 0.88 0.97 0.99 1.00     n=12  0.25 0.62 0.87 0.97 0.99 1.00 | 0.28 0.63 0.86 0.96 0.99 1.00     n=13  0.23 0.58 0.84 0.96 0.99 1.00 | 0.25 0.60 0.83 0.95 0.99 1.00     n=14  0.21 0.55 0.82 0.94 0.99 1.00 | 0.23 0.56 0.81 0.93 0.98 1.00     n=15  0.19 0.51 0.79 0.93 0.98 1.00 | 0.21 0.53 0.78 0.92 0.98 0.99     n=16  0.17 0.48 0.76 0.92 0.98 1.00 | 0.19 0.50 0.76 0.90 0.97 0.99     n=17  0.15 0.45 0.73 0.90 0.97 0.99 | 0.18 0.47 0.73 0.89 0.96 0.99     n=18  0.14 0.42 0.70 0.88 0.96 0.99 | 0.16 0.45 0.71 0.87 0.95 0.98     n=19  0.12 0.39 0.68 0.87 0.96 0.99 | 0.15 0.42 0.68 0.85 0.94 0.98     n=20  0.11 0.37 0.65 0.85 0.95 0.99 | 0.14 0.40 0.65 0.83 0.93 0.98     n=21  0.10 0.34 0.62 0.83 0.94 0.98 | 0.13 0.37 0.63 0.82 0.92 0.97     n=22  0.09 0.32 0.59 0.81 0.93 0.98 | 0.12 0.35 0.60 0.80 0.91 0.96     n=23  0.08 0.29 0.56 0.79 0.91 0.97 | 0.11 0.33 0.58 0.78 0.90 0.96     n=24  0.07 0.27 0.54 0.76 0.90 0.97 | 0.10 0.31 0.56 0.76 0.88 0.95     n=25  0.07 0.25 0.51 0.74 0.89 0.96 | 0.09 0.29 0.53 0.73 0.87 0.94     n=26  0.06 0.23 0.48 0.72 0.87 0.95 | 0.08 0.27 0.51 0.71 0.85 0.93     n=27  0.05 0.22 0.46 0.69 0.86 0.95 | 0.08 0.26 0.49 0.69 0.84 0.92     n=28  0.05 0.20 0.43 0.67 0.84 0.94 | 0.07 0.24 0.47 0.67 0.82 0.91     n=29  0.04 0.18 0.41 0.65 0.82 0.93 | 0.06 0.23 0.45 0.65 0.81 0.90     n=30  0.04 0.17 0.39 0.62 0.81 0.92 | 0.06 0.21 0.43 0.63 0.79 0.89


title: Determining if a criteria reduces scatter

text: I'm stuck on deciding what would be the best approach for my current problem: I have an array of approximately 1500 data points and I'm trying to find a best fit line through these points based on an complex expression with two free parameters. I found those parameters with Matlab's _lsqcurvefit_ function and I calculate the r squared to be approximately 0.7 after fitting the line. My next step is to impose a criteria to the original dataset, which reduces the number of valid data to ~ 400. Since the criteria is one of the basic assumptions for the theory behind the dataset, I expect to see a "reduction" in the scatter, i.e. increase in r squared. First thing that came to mind was doing the fitting procedure again on the reduced dataset and then comparing the two new fitting parameters with the old ones (theory doesn't say anything about the default or best values of the parameters). However, this doesn't say anything about the scatter. Regarding this, I've looked at this link, which describes using a t-test to compare the values of the parameters. However, the sample sizes need to be approximately the same which isn't true in my case (1500 <-> 400). Any advice or help is greatly appreciated.


title: Is a repeated measures ANOVA appropriate for multiple within-subject measurements?

text: I recorded 4 groups of captive birds (each group a difference species) at the same 3 times each day for 15 days. I want to determine if there is a difference in their vocalization properties such as call rate (the number of chirps per second for each call), call length, etc. during these 3 specific times of day. I repeated my recordings over several days because of the varying number of calls - sometimes the birds were silent, other times there were 10 calls in one recording session, etc.   1. My first reaction is to run a RM ANOVA (in SPSS), and if the differences between time of day are significant, do post-hoc tests for further investigation. Does this seem appropriate? I was initially going to do individual tests for each group, but now I am wondering if I can do the RM ANOVA using all of the groups and controlling for species differences.   2. I was thinking of just analyzing the actual calls that were recorded, ignoring the fact that some days there were no calls to analyze, as I am not interested in the differences between days, only between times. The problem with this is that I then have unbalanced data - for example, I may have a total of 22 calls to analyze at the beginning of the day, 17 at the middle, and 35 at the end. Would this be an issue with a RM ANOVA? I have been confusing and second-guessing myself like crazy - so I appreciate any advice!


title: Calculations using data captured over time seems to need a more advanced analysis method

text: I got a couple comments on my question. First of all, thank you very much. I was asked to be more clear. So, here's an attempt at that. I'm capturing raw data of play performance for each player and each game in the whole league. Then I can compute things like Goals Against averages and Shooting Percent for a selectable set of players over a selectable number of games, games going from the most recent and including the games that go back into the past for a selectable number. Then for one particular team I can compare the calculations against the same calculations for their opposing team, hoping to get an idea of what to expect by way of an outcome in an upcoming match between those two particular teams. The problem is that as the number of games increases the calculations for the two teams in question approach being equal to each other which starts to make it "look" like the teams are of the same quality when they're not in real life. When I said "noise" above I wasn't being clear... sorry. It's just that I wonder if a lot of what is really going on in real life IS captured in the data but the information is being obscured somewhere in all the data that is common between the two teams. For instance if the two teams have both played 40 games and one of the teams has won 38 and lost 2 and the other has won 35 and lost 5 the calculations over the entire 40 games will make the two teams "look" the same when really the 3 game difference in real life is very meaningful. I'm hoping someone could just kinda point me in a direction to study up on statistical techniques to formally and rigorously get the benefit of a long data sampling duration without making all the participants start to "look" the same. Certain kinds of correlation functions maybe? Time series techniques? Even just a little hint from some of you folks who have a good grip on this stuff would help me lots and lots!! Thanks for your comments!!!


title: How does reliability and validity affect the results (descriptive statistics)?

text: Okay, I think this might be a stupid question but I just can't quite seem to get there. I understand what reliability and validity is, but I don't understand just what I should do with it when looking at the results. Imagine a super-simple test examining how much people spend on groceries every week. For example, say that that test gives a significant result that is reliable but with low validity. So I get that that means that the test wasn't good for looking at the thing we wanted to look at, but that we would get roughly the same results if we repeated it. What can I do with that result then? Does it make the result completely useless? Anything else than draw the conclusion that we can't say anything about the subject that we wanted to look at? The opposite as well: a valid test with low reliability. Anything else than drawing the conclusion that we can't say anything? I'd be very grateful for any help. Thank you.


title: Expectation of a squared Gamma

text: If a Gamma distribution is parameterized with $\alpha$ and $\beta$, then: $$ E(\Gamma(\alpha, \beta)) = \frac{\alpha}{\beta} $$ I would like to calculate the expectation of a squared Gamma, that is: $$ E(\Gamma(\alpha, \beta)^2) = ? $$ I _think_ it is: $$ E(\Gamma(\alpha, \beta)^2) = \left(\frac{\alpha}{\beta}\right)^2 + \frac{\alpha}{\beta^2} $$ Does anyone know if this latter expression is correct? Thanks!


title: Probability Interval for Gamma Distribution

text: There is something that I am doing wrong in the exercise below and I would appreciate some help figuring it out. Let $X_1, X_2, \ldots, X_5$ be a random sample from a $\Gamma \left(3,3 \right)$ distribution . Define $W=\sum_{i=1}^5 X_i $.   1. Obtain the distribution of $W$ and of $2W/3$   2. Find $c_1$ and $c_2$ such that $P \left( c_1 < W <c_2 \right)=0.95 $ * * * Since $W$ is the sum of iid gamma variables, $W \sim \Gamma \left(15,3 \right)$ and now a straightforward transformation argument shows that $Y=\frac{2}{3} W \sim \chi^2 \left(30 \right) $. This completes question 1. Now for question 2, my thoughts were to transform these probabilities using the result from question 1, namely that $\frac{2}{3} W \sim \chi^2 \left(30 \right) $. Hence, unless I have made a horrible mistake here: $$P \left( \frac{2}{3} c_1 < Y < \frac{2}{3} c_2 \right)=0.95 $$ We can find the constants from a chisquare table for $30$ degrees of freedom. In my case, they are $16.791$ for the lower bound and $46.979$ for the upper one. Then if we equate $\frac{2}{3}c_1=16.791$ and solve for $c_1$ we should be able to find the constant for the above gamma confidence interval, right? The problem is that according to $R$, $P \left( W<c_1 \right) =1 $, instead of the desired $0.025$ which means I have gone wrong somewhere. I do not think it is in question though $1$ as I have verified the result many times. Thank you.


title: Analysis with binomial outcome and categorical predictor

text: I have a need to analyze an experiment that had a binomial outcome and categorical predictors. For context we sent out an email to our customers and encouraged them to register their account online (did register or did not register online). There are different nine groups that received the email. The groups are segmented based on their account opening date with our firm (1 week ago, 2 weeks ago, 3 weeks ago, etc...). What method would I use to determine if the response from each of the nine groups is significantly different? We have used R and ran ANOVA in the past when our dependent is continuous but I am not sure how to apply the same procedure with a binomial. Thank you in advance for any help!


title: Duplicate data for SVM

text: Can we use duplicate data as an input to SVM? The duplicate data that I mean is, let say we have 50 of same data (maybe being duplicate) from total of 100 data. Will this kind of data effect the performance result that we obtained from SVM? Do we have to avoid any duplicate data for SVM?


title: Calculating $Var\left\{(\hat{m}-m)^2\right\}$ for a univariate normal distribution

text: Suppose $\hat{m} = \frac{1}{N}\sum_{i=1}^{N}(X_i)$ where $X_i \sim N(m,\sigma)$. Are the following steps correct? $Var\left\{(\hat{m}-m)^2\right\} = E\left\{(\hat{m}-m)^4\right\} - E^2\left\{(\hat{m}-m)^2\right\}$ $= 3E^2\left\{(\hat{m}-m)^2\right\} - E^2\left\{(\hat{m}-m)^2\right\}$ $= 2E^2\left\{(\hat{m}-m)^2\right\}$ ~~and I know that $ E\left\{(\hat{m}-m)^2\right\} = \frac{1}{N^2}\sigma$.~~ (I was wrong here. Read the Update) Then, $Var\left\{(\hat{m}-m)^2\right\} = 2\frac{1}{N^4}\sigma^2$ * * * However the textbook says (without any proving) that $Var\left\{(\hat{m}-m)^2\right\} \tilde{} \frac{1}{N^2} $ Where am I going wrong? **Update:** as whuber told in the comments, i was wrong about $ E\left\{(\hat{m}-m)^2\right\} $. This expectation equla to $\frac{1}{N}\sigma$ and not $\frac{1}{N^2}\sigma$. Therefore, the variance is $Var\left\{(\hat{m}-m)^2\right\} = 2E^2\left\{(\hat{m}-m)^2\right\} = 2\frac{1}{N^2}\sigma^2 \tilde{} \frac{1}{N^2}$ Anyway, the answer provided by mpiktas is also correct and i prefer to chose it as the best answer.


title: Variation in PCA weights

text: I have weights of SNP variation (output through Eigenstrat program) for each SNP for the three main PCs. I wish to reduce my list of SNPs to those that show maximum differentiation between the three PCs. Can anyone help me with which statistical method to use to do this. say, if each PC describes the magnitude of variation, how can I find mutually exclusive rows or rows that define the strongest differenciator of each PC. e.g: > SNPNam PC1-wt PC 2-wt PC3-wt   >  SNP_1 -1.489 -0.029 -0.507   >  SNP_2 -1.446 -0.816 0.661   >  SNP_3 -1.416 0.338 1.631   >  SNP_4 -1.392 -1.452 0.062   >  SNP_5 -1.278 0.362 -1.006   >  SNP_6 -1.21 0.514 0.144   >  SNP_7 -1.119 -0.633 0.163   >  SNP_8 -1.112 -0.193 -0.256   >  SNP_9 -1.054 1.081 -1.519   >  SNP_10 -0.936 -1.052 -0.419   >  SNP_11 -0.861 0.381 -0.207   >  SNP_12 -0.662 0.852 -0.211   >  SNP_13 -0.503 -1.602 0.585   >  SNP_14 -0.417 0.529 -1.003   >  SNP_15 0.101 -0.85 -0.258   >  SNP_16 0.198 -0.435 -1.599   >  SNP_17 0.588 -1.292 -1.257   >  SNP_18 1.167 0.891 1.106   >  SNP_19 1.35 0.036 0.729   >  SNP_20 1.532 1.599 0.499 Any help regarding which test to perform and how to (which package) is very much appreciated.


title: How to draw a probable outcome from a distribution?

text: I have collected positional data. To visualize the data, I'd like to draw a 'typical' outcome of an experiment. The data comes from a few hundred experiments, where I identify a variable number of objects at different positions relative to the origin in 2D. Thus, I can calculate the average number of objects, as well as estimate the empirical distribution of the objects. A plot of the 'typical' outcome would then have the average (or possibly mode) number of objects, say, 5. What I'm not sure about is where to position these 5 objects. To simplify the problem, assume that the data follows a 2D normal distribution. If I were just to randomly draw 5 points from the distribution, I might get one point at [3,3], which would be a very rare outcome, and would thus not reflect the 'typical', or 'average' outcome. However, just drawing 5 points at [0,0] would also not make sense - even though [0,0] is the average position of the objects, 5 overlapping points are not an 'average' outcome of the process, either. In other words, how can I get a 'likely' draw from a distribution? * * * **EDIT** It looks like I should mention _why_ I don't want to use the usual methods (like a 2D smoothed histogram, or plotting all the many points) to look at the 2D distribution.   1. The objects (which are vesicles (i.e. little spheres) inside cells) vary in number, size and position (distribution of the distance from the cell center, amount of clustering). I would like to display all these features in one graph. Since there are several hundred cells containing many vesicles each, it is not very useful to combine them all in a single plot. I am well aware that I could use a multipanel graph showing the distributions of all parameters, but this would be a lot less intuitive.    2. I would like to show a 'typical' cell that shows all the salient features that characterize a specific phenotype. This way, if I want to image a particular phenotype in a mixed population, I know what kind of cell I'm looking for.   3. I think such a plot would be a cool way to display a lot of information at once, and I just want to try. Maybe it would be clearer If I said that I want to simulate a likely experimental result based on my measurements?


title: How to test for spatial correlation in the residuals of a dynamic panel data model?

text: I have an unbalanced panel on daily data for 8 months for 8 years on a pollutant. Thus, my T is very large and my N(number of cross-section) is 112. My model contains a lagged dependent variable as a regressor or what economists refer to as dynamic panel data model. In addition, I have year and time fixed effects and some other variables. I wish to test whether my residuals from an OLS regression are spatially correlated where the spatial correlation dies out as the distance between the cross-section units decreases. I don't know how can I do so in the dynamic panel framework with missing observations and a very large T.So please suggest something. Many thanks in advance.


title: Recommended method for finding archetypes or clusters

text: I wish to cluster users together in a database, with each user represented by a number of features that are both discrete and continuous in nature. The aim is to define a small number of archetypal "users" with specific set of features. All other users are then categorized as being similar to one or other of these archetypes. An important consideration is that I expect the features to have strong dependency structures, and I would like the method to be effective at making these explicitly visible. For example say the features per user are:   * gender (m/f)   * location (one of 10 cities)   * favorite color (red/green/blue). Let's say that we have N users and that favorite color is a R.V. dependent on gender and city. How are we to discover possible strong correlations with gender and/or location and favorite colors? There are a number of clustering techniques, from K-NN, k-means, matrix factorization, even PCA, but many seem to hide the underlying correlations that tie the users together. Could anyone recommend suitable methods for this unsupervised learning task? * * * [heavily edited in an effort to revive and resolve]


title: Fit lme in R with multiple levels of grouping and with correlation structure defined for more than inner-most level of grouping?

text: I’m trying to model the correlation/covariance between a pair of variables measured concurrently and repeatedly on a group of subjects. Since both variables are measured responses I want to avoid declaring one to be a response variable and the other a ‘fixed’ covariate. Here’s a synopsis of the data frame:               > names(SAmv)     [1] "Subj"     "State"    "ApoA1Ag"  "TlActVWF"          > table(SAmv$Subj)          13 25 27 30 34 40 46 47 48 51 55 57       6  4  3  3  5  3  3  4  3  5  6  3          > table(SAmv$State)          1  2  3 10 11 12 13 14 15       7  5  4  8  8  8  4  3  1          > dim(SAmv)     [1] 48  4      State is an investigator assigned factor which ‘measures’ the general state of the of the Subj at the time an observation is taken. I want to estimate the covariance between ApoA1Ag and TlActVWF at the Subj level and the the State within Subject level. To this end I have reshaped/reordered the data frame as rs2.SAmv (96 x 4, with variables: Subj, State, var, resp) where var is a label for which response variable is observed in resp and the data is ordered Subj/State/var . I have fit the model:               rs2.SAmv.lme2 <- lme(fixed = resp ~ var, data = rs2.SAmv ,     random = ~ var | Subj, weights = varIdent(form = ~1 | var),     correlation = corCompSymm(form = ~ 1 | Subj/State) )      This models the heteroscedasticity of the two response variables and gives an estimate of the correlation between ApoA1Ag and TlActVWF at the State within Subj level (inner most level of grouping). Do you have any suggestions for fitting a model which would also provide an estimate of the correlation between these two variables at the Subj level? It seems this will require a custom corStruct constructor function but I’m not conversant enough in R to write it.


title: How do I analyze a 2x2 experiment with DV measured pre & post intervention?

text: In my experiment, each subject is assigned to one of four conditions (2 IVs, 2 levels each). The DVs are measured pre and post intervention. Is this a 2x2 with repeated measures, or is it a 2x2x2, or something else? Perhaps it is a split plot design?


title: Which model to use with repeated measures data that contains multiple binary dependent variables

text: What model should I use??? I have daily repeated measures data. It has multiple dependent presence absence variables, (of which, I have collapsed into a CA with continuous variables of CA1 & CA2....or I can choose to not collapse the data and leave each dependent variable as a binary variable). My independent variables are weather data that was collapsed into a PCA (or I can use each individual weather measurement), 2 categorical variables, and time (or a covariate of time). I am wanting to create a model that allows me to solve for time (or the covariate of time). I have been looking at SEM latent growth models, logistic regression, time series, or mixed regression models. The software that is available to me is SAS, SPSS, PCord, and Canoco. "Gung":More information: CA: correspondence analysis. Presence and absence of insects on decomposing organisms that were measured periodically, mostly daily for the first 2 months and then periodically after that as insect presence greatly decreased. Two organisms were put out in the spring at different times, then some more in the summer, then fall, then winter. The time component that I am wanting to predict is how long that the organism has been decaying...so once I have a model of my data, I could put in all of the measurable variables of a decomposing organism that I find, and solve for time (decomposition day: how many days that the body has been decomposing). By covariate of time, which I think would actually be a better predictor is accumulated degree days (ADD). This is the sum of the daily min and max temp divided by two since the start of decomposition. This represents the accumulated heat units of the organism, and has been shown to influence different species of insects. EX) typically, flies are the first to a decomposing organism, so at smaller ADD values. Later, beetles appear at large ADD values. All of my independent variables have had literature published that suggest that they influence insect occurrence. As a second part, I would like to look at the rate of change in community structure (or insect occurrence) over time...or meaning as decomposition progresses. Thank you so much for any advice!!


title: Suitable data preprocessing for an image recognition task: normalization ,standardization or neither?

text: I aim to train my classifier for an image recognition task. What kind of preprocessing steps I need to take to enhance my results? (I have >40000 images with >700 pixels so large amounts of data with a vast feature vector.) Also can you please suggest any other efficient preprocessing method to advance performance after you share about normalization, standardization?


title: Using non-treatment factors to rank patients for receiving a treatment

text: I'm facing a problem and am wondering if the biostats community may already have a method for it that I just haven't heard about in my field. I am doing some work with a test vs. control group and a certain treatment. Measuring the treatment effect by itself (between the two groups, overall) I am fine with. But the ability to give the treatment is scarce, and in the future when we select certain individuals out of a pool to administer to, we want to only give it to those with the largest expected change of a certain variable we are interested in (call it Y). So yes, on average we know how much of a change to expect thanks to the treatment, but this doesn't help learn which demographic feature(s) we should use to choose future individual patients by, if our goal is to maximize each individual's change in Y (relative to a similar patient in control; if they were going to have that change in Y without treatment we certainly don't want to use it on them). So I'd like to determine, after the basic fact of yes the treatment was given, which other demographic factors are associated with a change in Y that is larger than what would be found in the control group. For example, I'd like to be able to say some sentence like, "oh patients with high values on variable X4 and X7, in addition to getting treatment, have a much greater change in Y than we would expect to see in a control group with similar demographic values". Is there some known method in biostats that exists for this? Many thanks for your help.


title: Relationship between threshold, number of features and accuracy/error? (text classification)

text: Are you aware of any research papers that explain a relationship between the following concepts?   * threshold (removal of features, whose frequencies are greater than or less than a defined threshold)   * number of features   * error


title: Non-normality in MANOVA with significant planned contrasts versus non-significance in Kruskal-Wallis- A Conundrum

text: I am conducting a between groups analysis for 3 groups ( _N_ =26, one control group and two experimental groups) with 4 dependent variables. I initially planned to conduct a MANOVA, however I am not meeting normality (significant skewness and kurtosis using _z_ > 2.58) and homogeneity of variance assumptions for 2 of the DVs. Following a log transformation, one DV was corrected by the other has remained significantly skewed. I initially decided to conduct a Kruskal-Wallis non-parametric test. However, due to a non-significant main effect I was unable to carry out post hoc contrasts. The initial MANOVA I conducted was significant on all DVs for several planned contrasts (control vs both experimental groups- one of my study's main hypotheses), however it seems I will be unable to report these findings if I choose a non-parametric test as above. Is there any way I can   * Report these contrasts in MANOVA despite non-normality?   * Conduct planned contrasts in Kruskal-Wallis (which does not seem like an option according to Field, 2004)   * Use an alternative analysis?


title: Tutorial for performing Contrasts for MANOVA

text: ### Context I have a multivariate dataset with a test group and three control groups. I was thinking that the best way to determine if and how the test group differed from all of the control groups would be to perform a MANOVA, then perform contrast analysis between the test group and all of the control groups. ### Main Question I cannot seem to find a good description of how contrast analysis for MANOVA is performed.   Every tutorial/online class I can find, just talks about throwing the contrast into SPSS, or the univariate case.   * Can anyone point me to a good tutorial on contrasts for MANOVA?  ### Questions   * Is MANOVA with contrasts the appropriate way of answering my research question?   * Would it be better to perform repeated Hotelling's tests on the test group vs each control group, and just use a Bonferroni correction?


title: How to extract data from published articles (RCTs) to do a meta-analysis?

text: I am doing a meta analysis for the first time and have a few basic questions regarding statistical analysis. Let's say I have one study where the primary outcome (thrombosis) in the 2 treatment groups (intervention v. placebo) was compared using the Mantel-Haenszel chi-square test. It does not report df in the article (is this something I can deduce?). Then it tells me that "Treatment effects were expressed as a weighted average of the strata-specific relative risks." Is this something I need to consider for the meta-analysis statistical calculation? If so, what does it mean? In the results we learn "Among the 866 participants who had patency assessed, the primary outcome of fistula thrombosis at 6 weeks occurred in 53 participants (12.2%) in the clopidogrel group compared with 84 participants (19.5%) in the placebo group (relative risk, 0.63; 95% CI, 0.46-0.97;P=.018). What do I do with all this information? I tried using a computerized program, but I didn't know what to fill all the blanks in with... In the second study, we have the following: "The hazard ratio was 0.81 (95% CI, 0.47 to 1.40) in favor of aspirin and clopidogrel therapy, but the reduction was not significant (P= 0.45). Although the event rates in the two treatment arms converged, there was no evidence that the proportional-hazards assumption was violated (P 0.53). The annual hazard rate for thrombosis among articipants receiving placebos was 0.59 (95% CI, 0.39 to 0.87), compared with 0.47 (95% CI, 0.31 to 0.71) among participants receiving aspirin and clopidogrel. The absolute risk reduction for the group receiving aspirin and clopidogrel was 0.03 (95% CI,-0.19 to 0.07), and the number needed to treat was 33.3 (95% CI, 14.9 to >/=1000)." Regarding analysis, it states: "The cumulative incidence of the first episode of thrombosis was estimated with the Kaplan-Meier method, and differences in rates between treatment groups were tested with the log rank test. The cumulative incidence of the first bleeding event was analyzed in a similar manner." What do I do with all this information? Thank you for any help you can provide.


title: Monte carlo integration in spherical coordinates

text: I was playing around with writing a code for Montecarlo integration of a function defined in spherical coordinates. As a first simple rapid test I decided to write a test code to obtain the solid angle under an angle $\theta_m$. For two random number $u$ and $v$ in $[0,1)$. I generate a homogeneous random sampling of the spherical angle using $$\phi=2\pi u$$ $$\theta = \arccos(1-2v)$$ For $N$ generated points, I have $M$ points for which $\theta < \theta_m$. My first idea was that since I have an homogeneous sampling I should have obtained the correct solid angle $\Omega=2 \pi (1-\cos (\theta_m))$ simply as $4\pi\times M/N$. Actually it looks like I the correct result comes out only if I use: $$\Omega=\sum_{i=1}^M \frac{4\pi}{N} 2 \cos(\theta_i)$$ I can not see the reason why this should be correct. The probability distribution function in $\theta$ is $PDF=1/2 \sin(\theta)$ so I would rather expect I should normalize each point of the sum by this function but this doesn't work. What am I doing wrong and how could I justify the cosine? Many thanks!


title: Sample size formula for an F-test?

text: I am wondering if there is a sample size formula like Lehr's formula that applies to an F-test? Lehr's formula for t-tests is $n = 16 / \Delta^2$, where $\Delta$ is the effect size ( _e.g._ $\Delta = (\mu_1 - \mu_2) / \sigma$). This can be generalized to $n = c / \Delta^2$ where $c$ is a constant that depends on the type I rate, the desired power, and whether one is performing a one-sided or two sided test. I am looking for a similar formula for an F-test. My test statistic is distributed, under the alternative, as a non-central F with $k,n$ degrees of freedom and non-centrality parameter $n \lambda$, where $\lambda$ depends only on population parameters, which are unknown but posited to take some value. The parameter $k$ is fixed by the experiment, and $n$ is the sample size. Ideally I am looking for a (preferably well-known) formula of the form $$n = \frac{c}{g(k,\lambda)}$$ where $c$ depends only on the type I rate and the power. The sample size should satisfy $$ F(F^{-1}(1-\alpha;k,n,0);k,n,n\lambda) = \beta,$$ where $F(x;k,n,\delta)$ is the CDF of a non-central F with $k,n$ d.o.f. and non-centrality parameter $\delta$, and $\alpha, \beta$ are the type I and type II rates. We can assume $k \ll n$, _i.e._ $n$ need be 'sufficiently large.' My attempts at fiddling with this in R have not been fruitful. I have seen $g(k,\lambda) = \lambda / \sqrt{k+1}$ suggested but the fits have not looked very good. **edit:** originally I had vaguely stated that the non-centrality parameter 'depends' on the sample size. On second thought, I found that too confusing, so made the relationship clear. Also, I can compute the value of $n$ exactly by solving the implicit equation via a root finder ( _e.g._ Brent's method). I am looking for an equation to guide my intuition and for use as a rule of thumb.


title: Getting the weights to compare AMORE NN models?

text: I'm using the AMORE NN package to build multiple models from the same training set and i want to compare the weights between the models. I'm pretty much an R newbie, but some of the stronger R users I know can't answer either. here's the code:               library("AMORE")     net <- newff(n.neurons=c(104,200,104), learning.rate.global=1e-2,                   momentum.global=0.5,error.criterium="LMS", Stao=NA,                      hidden.layer="tansig",output.layer="purelin",                   method="ADAPTgdwm")     result <- train(net, x, y, error.criterium="LMS",                      report=TRUE, show.step=50, n.shows=20 )      `net` is pretty easy to break out into some functions and I think input weights. `net$neurons` is where all the layer weights are. `net$neurons` is a or list of length 304. `net$neurons[n]` will print out what appears to be a set of units from the model. But it contains many things and I can't seem to get them out short of cut and paste. I won't paste it all here - my model is pretty big as you can see and it goes on for several pages, but just to give you an idea: > net$neurons[1] > > $type > > [1] "output" > > $activation.function > > [1] 3 > > $output.links > > [1] NA > > $output.aims > > [1] 103 > > $input.links [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 [19] 19 20 21 > 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 [37] 37 38 39 40 41 42 43 44 45 > 46 47 48 49 50 51 52 53 54 [55] 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 > 70 71 72 [73] 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 [91] 91 > 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 [109] 109 110 > 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 [127] 127 > 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 [145] > 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 > [163] 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 > 180 [181] 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 > 197 198 [199] 199 200 > > $weights [1] 7.767180e-03 -5.118608e-03 4.683484e-03 6.499882e-03 > 6.897404e-04 [6] 4.169454e-03 -5.282092e-05 9.589098e-04 7.031523e-03 > 6.669872e-03 [11] 2.716605e-03 -4.188492e-04 -4.141848e-03 -6.576378e-03 > -6.812882e-03 [16] -3.511420e-03 -3.091805e-04 7.128892e-03 7.160105e-03 > -8.050284e-03 [21] 3.386642e-03 2.579751e-03 7.065484e-03 1.872059e-04 > 1.829110e-03 [26] 7.472301e-03 -8.303545e-03 7.981113e-03 3.584151e-03 > -4.432509e-03 [31] -8.705997e-04 -7.456132e-03 -9.884449e-04 7.523373e-03 > -6.673688e-03 [36] -9.070783e-04 3.068912e-03 -8.045977e-04 -2.165867e-03 > 5.507719e-03 ... There's lots more, lots more, but it all has this general format.   Things that give no answer at all:               net$neuron[1]$weights     net$neuron[1,"weights"]     net$neuron[1,1]     net$neuron[1]["weights"]      probably some of the above is laughably impossible anyway, but can anyone point me somewhere useful? thanks


title: Using hierarchical clustering to classify?

text: Given an hierarchical clustering of data points, some of which are labeled, are there good ways to use the tree/dendrogram to make predictions for the unlabeled points? One approach might be to find the "best" place to cut the tree so that clusters match labels. I'd be especially interested in efficient ways to cut each cluster at a different height. But I'd also be interested in approaches that don't use "hard" cuts to make the predictions.


title: Problem with zero epsilon value for Kalman filter

text: I am coding a Kalman filter in R and I am trying to estimate the parameters dV and dW using the methods as provided in Vignette section2. However I get Epsilon=0 (which is dV for Kalman filter) for my ts data. I am not sure how we get a non-zero dV for the Nile data (available in R) used in the example in 1. Can somebody explain what dV means as it is not the same as var(z)?               > dput(z)   # my time series data     structure(c(4341, 3426, 2312, 2118, 2003, 1881, 1734, 1764, 1634,      1649, 1316, 1191, 1145, 1094, 1118, 1107, 1145, 1087, 1071, 1056,      1075, 984, 925, 963, 916, 873, 836, 837, 903, 924, 847, 863,      874, 872, 926, 1017, 1274, 1306, 1334, 1240, 1971, 2620, 2803,      2701, 2846, 2728, 2597, 2440, 2673, 2588, 2578, 2899, 3129, 3100,      4244, 5175, 6409, 7444, 7732, 5847, 4483, 3350, 2907, 2626, 2162,      2151, 1628, 1575, 1353, 1443, 1150, 1161, 1122, 1147, 1173, 1026,      966, 925, 980, 1002, 939, 1010, 958, 1141, 1138, 1112, 1029,      1369, 1878, 1478, 1694, 2065, 2635, 2438, 2539, 3099, 3014, 2825,      2610, 2972, 2911, 3508, 4579, 5421, 4382, 3997, 3927, 4116, 4122,      3696, 3793, 3607, 3386, 2817, 2293, 2101, 2111, 1933, 1769, 1694,      1711, 1560, 1428, 1322, 1286, 1346, 1174, 1006, 1031, 997, 1066,      1071, 1015, 1064, 1004, 1003, 992, 1095, 1161, 1252, 1416, 1495,      1649, 1918, 1913, 1959, 2341, 2609, 2952, 2918, 2684, 2836, 2873,      3205, 3078, 3510, 3524, 3126, 2685, 3213, 4183, 5710, 5925, 6565,      5543, 4649, 3324, 2673, 2156, 1899, 2175, 1892, 1693, 1539, 1411,      1381, 1339, 1385, 1359, 1336, 1310, 1178, 1204, 1186, 1224, 1200,      1209, 1308, 1212, 1172, 1323, 1495, 1580, 1758, 1743, 2019, 2267,      2363, 2522, 2680, 2889, 2894, 2386, 2934, 2908, 2615, 2434, 3022,      4408, 4546, 4965, 7135, 9531, 9758, 8440, 7945, 6356, 4819, 4100,      3085, 2694, 2202, 2023, 1828, 1812, 1518, 1404, 1427, 1361, 1221,      1174, 1070, 1059, 969, 982, 1059, 1009, 1011, 999, 969, 1042,      1004, 1133, 1249, 1380, 1519, 1746, 1907, 2053, 2173, 2412, 2764,      2742, 2990, 3085, 2576, 3108, 2990, 2972, 3360, 3625, 3257, 3375,      4520, 5707, 8206, 7994, 6058, 5132, 4153, 3054, 2427, 2212, 2050,      1655, 1755, 2048, 4878, 3341, 2561, 2464, 2158, 1881, 1843, 1864,      1970, 1799, 1515, 1473, 1793, 1878, 1919, 1854, 1977, 2500, 3608,      5708, 8518, 10915, 13458, 12567, 10789, 8682, 6933, 6048, 4972,      4132, 3323, 3362, 3145, 2976, 3181), .Tsp = c(2004, 2009.98076923077,      52), class = "ts")          > StructTS(z,"level")           Call:     StructTS(x = z, type = "level")          Variances:       level  epsilon        403733        0            >data(Nile)     >StructTS(Nile,"level")          Call:     StructTS(x = Nile, type = "level")          Variances:       level  epsilon          1469    15099


title: What is $\text{cov}$$(\hat{y_i}, \hat{y_i}^*)$?

text: Using the simple linear regression model: $ {y_i}= {\beta_0} + {\beta_1}x_i + \epsilon_i$, where E[$\epsilon_i$]=0 and var[$\epsilon_i] = \sigma^2$... If $ \hat{y_i}= \hat{\beta_0} + \hat{\beta_1}x_i$ and $ \hat{y_i^*}= \hat{\beta_0} + \hat{\beta_1}x_i^*$ (where $\hat{\beta_0}$ and $\hat{\beta_1}$ are ordinary least squares estimates of $\beta_0$ and $\beta_1$) how would I go about deriving $\text{cov}$$(\hat{y_i}, \hat{y_i}^*)$? Here is what I have been thinking: $\text{cov}$$(\hat{y_i}, \hat{y_i}^*) = cov(\hat{\beta_0} + \hat{\beta_1}x_i, \hat{\beta_0} + \hat{\beta_1}x_i^*) = \beta_1^2 cov(x_i, x_i^*) = \beta_1^2[E[(x_i - E(x_i))(x_i^* - E(x_i^*))]] $ If the work that I've done is correct, is there any more steps that I can take?


title: Meaning of risk

text: Suppose you want to look at the risk of glaucoma in 7 groups (e.g. grouped by weight: 100-110 lbs, 111-120 lbs, etc..). Is risk measured by the proportion of people that have glaucoma? Or should you look at excess risk/relative risk? I guess it shouldn't matter what you look at because the risk pattern would be the same?


title: Degrees of freedom for Kuiper's Test of uniformity

text: How do you determine the degrees of freedom for Kuiper's Test of Uniformity for circular data? I have been using various R packages for circular data (circular, CircStats) and I can't seem to find out how to get the degrees of freedom for the Kuiper's test.


title: Does it make sense to calculate Q2 and R2 values on PLS-DA models?

text: Since PLS-DA is a computational technique which deals with outcomes expressed as a categorical variable (e.g. "Yellow","Brown","Black","Green") I cannot understand how it is possible to calculate Q2 and R2 values. In facts, these values are calculated based on the difference between the real outcome and the outcome predicted by the model after removal of some features; this procedure seems correct with continuous outcomes but how is it possible to perform a difference between "Green" and "Yellow" for instance? Many thanks.


title: Prediction using CRFs for time series data

text: I have a little confusion about validity of some predictions I am making using a CRF model I have trained. The CRF model is trained on some input time-series, and when making predictions, I am passing the entire sequence to be labelled to the model. I have a concern that this is incorrect for the problem, since what I want is predictions for each observation in the time series up to that point (as if it were coming in, in real time), and I have the feeling that predictions are being made using the entire sequence (mainly due to references I see about a forward-backward algorithm for inference). Do I need to pass only the sequence up to "what has currently been observed" in order to get a causally predicted label for that point in time, or is there such a thing as "forward only" prediction for such scenarios?


title: What is the correct notation if two random variables belong to the same distribution?

text: I want to explain that both the real and imaginary part of a complex variable follow a zero-mean complex Gaussian distribution. How can I write that? For one variable I'd write $a \sim \mathcal{N}(0,\sigma^2)$. But for two? Maybe $\{a,b\} \sim \mathcal{N}(0,\sigma^2)$? What's the correct notation?


title: What is the term for a time series regression having more than one predictor?

text: It's pretty tough to search the Web for info on something when you don't know what words are commonly used to describe it. In this case, I'm wondering what it's called when you include another predictor in a time series. As an example, say I'm modeling a variable $X$ using AR(3): $ X_t = \varphi_1 X_{t-1} + \varphi_2 X_{t-2} + \varphi_3 X_{t-3} + \varepsilon_t $ I want my model to include the effects of another variable—say, $Y$—so my model is now described as: $ X_t = \varphi_1 X_{t-1} + \varphi_2 X_{t-2} + \varphi_3 X_{t-3} + \beta_1 Y_{t-1} + \beta_2 Y_{t-2} + \beta_3 Y_{t-3} + \varepsilon_t $ What would the term (or terms) be that distinguishes the former model from the latter?


title: create graphic with multiple data-columns in different appearance in jmp

text: in jmp 11 on a mac I have data I want to plot nicely in a 2d-diagram. For the datapoints I made a regression to get a function for a smooth curve, now I want to plot the datapoints and the function in one graph. However in the graphic-creating-window I can only select the appearance for all columns, i.e. either 'dots' (for the datapoints) or 'function' (for the smooth function). If I select 'dots', both datacolumns are represented as only dots, if I select 'function', only the smooth function is shown. But I want to see the datapoints as dots and the function as continous line. Is this possible in jmp? Regards, Paul


title: Number of significant figures to put in a table?

text: Is there a well founded rule for the number of significant figures to publish? Here are some specific examples / questions:   * Is there any way to relate the number of significant figures to the coefficient of variation? For example, if the estimate is 12.3 and the CV is 50%, does that mean that the information represented by '.3' approaches zero?    * If a confidence interval has a range of orders of magnitude, should they still have the same number of significant figures, e.g.: 12.3 (1.2, 123.4) vs 12 (1.2, 120)   * Should the number of significant figures in an error estimate be the same or less than the number of significant figures in a mean?


title: How to model time-series temperature data at multiple sites as a function of data at one site?

text: I am new to time series analysis, and would appreciate any suggestions on how best to approach the following time-series regression problem: I have hourly temperature measurements at approximately 20 locations across one site over three years, along with static ancillary information (slope, elevation, aspect, canopy cover). The site is several hectares in size, and the temperature recording devices are spread across the site along a couple of transects, at ~20-50 m intervals. About 1 km away, I have hourly data from a weather station, which also provides measurements of wind speed, wind direction, humidity, solar illumination, etc. I would like to be able to predict the temperature (min,max,mean) at the site (in general) using only the data from the weatherstation; it is in place semi- permanently, whereas the temperature recorders at the site were only in place for 3 years. So in essence I have multiple independent variables (temperature, humidity, wind, etc) at one location (the weatherstation), but a single dependent variable (temperature) at multiple locations, each of which also has several time-invariant attributes: slope, elevation, aspect, etc. I am most interested in predicting the daily lows and highs at the site in general, rather than hourly temperatures at each temperature recording location in the site. Although, those hourly predictions would certainly be of value. My initial approach has been to compute daily average, minimum, and maximums from the temperatures at the site, and use these as dependent variables in simple linear regressions, using the measurements available at the weatherstation as independent variables. This works reasonably well (R2 > 0.50 with 2 predictors), but seems rather too simplistic for many reasons, and I imagine there must be more sophisticated (and powerful) ways to do this. For one, I'm not doing anything explicit about the time-series nature of the daily values in the regression, and although the min or average temp from one day to the next may not be as correlated as it is from one hour to the next, I wonder about issues with the independence of these daily data (or certainly hourly, if I were trying to predict hourly temperatures). Second, due to concerns with having multiple somewhat-correlated temperature measurements across the site (they are much more similar among themselves than any are to the weather station data), I am simply using the mean or min or max of all measurements across the site, versus including the data from each individual measurement location directly. But this also prevents me from using the time- invariant ancillary information from each temperature measurement location (slope, elevation, aspect, canopy cover), which presumably will explain a good part of the differences in temperatures between locations at the site. Third, due to concerns with the regression being dominated by the very strong diurnal cycle in temperatures, I'm only looking at daily values instead of hourly. Any suggestions on better ways to go about this (especially in R), or where to start looking, would be most appreciated! I realize there are alot of R packages that deal with time-series, but I'm having trouble finding the best place to start with this type of problem as none of the examples I've seen really seem to reflect the situation I'm trying to model here. Update: thinking about this a bit more, it is not clear to me whether time- series models are really appropriate here because I am not interested in predicting what will happen at some future specific point in time. Rather, I'm simply interested in how temperatures at the site are related to temperatures (and other environmental variables) at the weatherstation. I thought that perhaps time-series analysis would be of value because I was concerned that subsequent temperature measurements might not be sufficiently independent. Certainly, one hour's temperature depends a great deal on the previous hour, but the dependence is weaker for daily data. In either case, is the time- correlation/non-independence of time-series data a valid concern that should be addressed if one is not interested in a time-series prediction?


title: Repeated measures with correlated measures (not time)

text: Most repeated measures ANOVAs have time as the repeated measure; I was wondering about using a repeated measure that is not time. Say we fed two groups of animals different diets. At the end of the experiment, we sample the tissues, and measure ~30 different compounds (e.g. different fatty acids [FA]). Animals are sampled but once. Each FA is not independent of each other, as some of these compounds are converted from one form to another. As such, they are frequently moderately correlated with each other. Would it be fair to treat each compound as a repeated measure, with diet as a between subjects factor? Thus, the interaction between Diet X FA would tell me if the FA content differed among diets? I note that in many papers, researchers would perform 25-30 separate ANOVAs; one on each compound. Yet, these compounds are not independent of each other, as they are each measured simultaneously on the same animals. Thanks for any pointers.


title: Optimal Rotation Before Curve Fitting

text: I am fitting a 1-dimensional polynomial to a set of data points $\{(x_1,y_1,z_1),...,(x_N,y_N,z_N)\} \subset \mathbb{R}^3$ using least squares regression. My first instinct was to choose $x$ as an explanatory variable and perform regression on the sets $\{(x_1,y_1),...(x_N,y_N)\}$ and $\{(x_1,z_1),...,(x_N,z_N)$ separately. This would yield two functions $f_y, f_z$. My curve would then be given by $\{(x,f_y(x),f_z(x)) : x \in [a, b]\}$. Where $a,b$ are appropriately chosen. However, I don't think this is optimal as the x-coordinates may be identical in a worst-case scenario. My instinct is that I should perform PCA on the data-set beforehand and use the principal axis as my explanatory variable. Is this optimal in some sense or are there better ways to do this? Thanks in advance.


title: Need to find appropriate statistical test for a longitudinal intervention study

text: The study: Comparing two groups (intervention vs control) across time using a behavioral measure as the outcome variable of interest. I will have baseline behavioral measurements as well as more than 2 time points throughout treatment. I'm interested in finding any significant improvements from baseline across time (within each group) and whether or not there were significant differences between groups. Would this be a repeated measures ANOVA or MLM? Would I need to control for baseline behavioral measurement? Thank you in advance for any help with this!!


title: What interpretation do the parameters of a generalised linear model with effect coding have?

text: library(lme4)         out <- glmer(cbind(incidence, size - incidence)                      ~ period                      + (1 | herd),                      data = cbpp,                      family = binomial,                      contrasts = list(period = "contr.sum"))          summary(out)     Fixed effects:                 Estimate Std. Error z value Pr(>|z|)     (Intercept) -2.32337    0.22129 -10.499  < 2e-16 ***     period1      0.92498    0.18330   5.046 4.51e-07 ***     period2     -0.06698    0.22845  -0.293    0.769     period3     -0.20326    0.24193  -0.840    0.401      I was never in a situation where I needed to fit a generalised linear model with effect coding (`contr.sum` for `R` users). Can I apply the same interpretation as in the linear model case? In a normal linear model the intercept would be the grand mean and the $\beta$s (parameters for `period1`, `period2`, `period3` and `period4 = (Intercept) - period1 - period2 - period3` the effects i.e. how the factor levels deviate from the grand mean. Here is how I think the analogous interpretation for generalised linear models goes. (I will exponentiate all parameters and hence transform the log- odds(-ratios) to odds(-ratios).) The intercept $\exp((\text{Intercept}))$ would then be the overall **odds** of success vs. failure (sticking here to classical binomial terminology) and the $\beta$s the **log-odds-ratios**. And we get the **odds** for e.g. `period1` by adding $\text{(Intercept)}+\text{period1}$ and then exponentiating: $\exp(\text{(Intercept)}+\text{period1})$. Is the $\text{(Intercept)}$ really the overall/medium **odds** and the $\beta$s **odds-ratios**?


title: Limiting Distribution of a Polya-Eggenberger Urn with Stopping Rule

text: I am trying to find literature on the limiting distribution of a Polya- Eggenberger Urn with a stopping rule contingent on the proportion of the two colors in the urn. To clarify, let's say that the urn has white and black balls and whoever is pulling balls out of the urn stops once the ratio of the balls in the urn is say $0.5$. My guess is that in the limit the ratio of white to black balls will always be $>0.5$. I have however no idea how to prove this and I expect that this has been proven in the past. I could simulate, but what I am looking for is a proof. Best, Dave


title: Scheduling Algorithm for a multi-server queue problem

text: I have 4 servers, n customers and m reports. At any time, a customer may request one of m reports. There are only 4 servers which are capable of generating reports. Each server can only process one report at a time. Over time I can record the average generation time for each customer which is independent of the server. I assume that the arrival is Poisson Process. I would like to come up with the optimal schedule minimizing the average queue size. In other words, if I currently have say 10 customers in the system requesting 10 different report, how do I pick the order of serving these customers (other than first come, first serve) I'm familiar with some Stochastic Processes and Probability Methods. However, I'm struggling with:   1. How do I determine the distribution if it depends on each report?   2. If I don't assume Poisson arrival, how do I determine the distribution?   3. What kind of algorithms would be the best to minimize the average time?


title: Excluding the scatter points from a feature

text: I have a set of data points that are supposed to sit on a locus and follow a pattern, but there are some scatter points from the main locus that cause uncertainty in my final analysis. I would like to obtain a neat locus to apply it later for my analysis. The blue points are more and less the scatter points that I want to find and exclude them by a sophisticated way without doing it manually.![enter image description here](http://i.stack.imgur.com/GHn3o.png) I was thinking of using something like `Nearest Neighbors Regression` but I am not sure whether it is the best approach or I am not very familiar how it should be implemented in order to give me an appropriate result. By the way, I want to do it without any fitting procedure. The transposed version of the data is a following:               X=array([[ 0.87 , -0.01 ,  0.575,  1.212,  0.382,  0.418, -0.01 ,  0.474,              0.432,  0.702,  0.574,  0.45 ,  0.334,  0.565,  0.414,  0.873,              0.381,  1.103,  0.848,  0.503,  0.27 ,  0.416,  0.939,  1.211,              1.106,  0.321,  0.709,  0.744,  0.309,  0.247,  0.47 , -0.107,              0.925,  1.127,  0.833,  0.963,  0.385,  0.572,  0.437,  0.577,              0.461,  0.474,  1.046,  0.892,  0.313,  1.009,  1.048,  0.349,              1.189,  0.302,  0.278,  0.629,  0.36 ,  1.188,  0.273,  0.191,             -0.068,  0.95 ,  1.044,  0.776,  0.726,  1.035,  0.817,  0.55 ,              0.387,  0.476,  0.473,  0.863,  0.252,  0.664,  0.365,  0.244,              0.238,  1.203,  0.339,  0.528,  0.326,  0.347,  0.385,  1.139,              0.748,  0.879,  0.324,  0.265,  0.328,  0.815,  0.38 ,  0.884,              0.571,  0.416,  0.485,  0.683,  0.496,  0.488,  1.204,  1.18 ,              0.465,  0.34 ,  0.335,  0.447,  0.28 ,  1.02 ,  0.519,  0.335,              1.037,  1.126,  0.323,  0.452,  0.201,  0.321,  0.285,  0.587,              0.292,  0.228,  0.303,  0.844,  0.229,  1.077,  0.864,  0.515,              0.071,  0.346,  0.255,  0.88 ,  0.24 ,  0.533,  0.725,  0.339,              0.546,  0.841,  0.43 ,  0.568,  0.311,  0.401,  0.212,  0.691,              0.565,  0.292,  0.295,  0.587,  0.545,  0.817,  0.324,  0.456,              0.267,  0.226,  0.262,  0.338,  1.124,  0.373,  0.814,  1.241,              0.661,  0.229,  0.416,  1.103,  0.226,  1.168,  0.616,  0.593,              0.803,  1.124,  0.06 ,  0.573,  0.664,  0.882,  0.286,  0.139,              1.095,  1.112,  1.167,  0.589,  0.3  ,  0.578,  0.727,  0.252,              0.174,  0.317,  0.427,  1.184,  0.397,  0.43 ,  0.229,  0.261,              0.632,  0.938,  0.576,  0.37 ,  0.497,  0.54 ,  0.306,  0.315,              0.335,  0.24 ,  0.344,  0.93 ,  0.134,  0.4  ,  0.223,  1.224,              1.187,  1.031,  0.25 ,  0.53 , -0.147,  0.087,  0.374,  0.496,              0.441,  0.884,  0.971,  0.749,  0.432,  0.582,  0.198,  0.615,              1.146,  0.475,  0.595,  0.304,  0.416,  0.645,  0.281,  0.576,              1.139,  0.316,  0.892,  0.648,  0.826,  0.299,  0.381,  0.926,              0.606],            [-0.154, -0.392, -0.262,  0.214, -0.403, -0.363, -0.461, -0.326,             -0.349, -0.21 , -0.286, -0.358, -0.436, -0.297, -0.394, -0.166,             -0.389,  0.029, -0.124, -0.335, -0.419, -0.373, -0.121,  0.358,              0.042, -0.408, -0.189, -0.213, -0.418, -0.479, -0.303, -0.645,             -0.153,  0.098, -0.171, -0.066, -0.368, -0.273, -0.329, -0.295,             -0.362, -0.305, -0.052, -0.171, -0.406, -0.102,  0.011, -0.375,              0.126, -0.411, -0.42 , -0.27 , -0.407,  0.144, -0.419, -0.465,             -0.036, -0.099,  0.007, -0.167, -0.205, -0.011, -0.151, -0.267,             -0.368, -0.342, -0.299, -0.143, -0.42 , -0.232, -0.368, -0.417,             -0.432,  0.171, -0.388, -0.319, -0.407, -0.379, -0.353,  0.043,             -0.211, -0.14 , -0.373, -0.431, -0.383, -0.142, -0.345, -0.144,             -0.302, -0.38 , -0.337, -0.2  , -0.321, -0.269,  0.406,  0.223,             -0.322, -0.395, -0.379, -0.324, -0.424,  0.01 , -0.298, -0.386,              0.018,  0.157, -0.384, -0.327, -0.442, -0.388, -0.387, -0.272,             -0.397, -0.415, -0.388, -0.106, -0.504,  0.034, -0.153, -0.32 ,             -0.271, -0.417, -0.417, -0.136, -0.447, -0.279, -0.225, -0.372,             -0.316, -0.161, -0.331, -0.261, -0.409, -0.338, -0.437, -0.242,             -0.328, -0.403, -0.433, -0.274, -0.331, -0.163, -0.361, -0.298,             -0.392, -0.447, -0.429, -0.388,  0.11 , -0.348, -0.174,  0.244,             -0.182, -0.424, -0.319,  0.088, -0.547,  0.189, -0.216, -0.228,             -0.17 ,  0.125, -0.073, -0.266, -0.234, -0.108, -0.395, -0.395,              0.131,  0.074,  0.514, -0.235, -0.389, -0.288, -0.22 , -0.416,             -0.777, -0.358, -0.31 ,  0.817, -0.363, -0.328, -0.424, -0.416,             -0.248, -0.093, -0.28 , -0.357, -0.348, -0.298, -0.384, -0.394,             -0.362, -0.415, -0.349, -0.08 , -0.572, -0.07 , -0.423,  0.359,              0.4  ,  0.099, -0.426, -0.252, -0.697, -0.508, -0.348, -0.254,             -0.307, -0.116, -0.029, -0.201, -0.302, -0.25 , -0.44 , -0.233,              0.274, -0.295, -0.223, -0.398, -0.298, -0.209, -0.389, -0.247,              0.225, -0.395, -0.124, -0.237, -0.104, -0.361, -0.335, -0.083,             -0.254]])


title: How to understand the Gaussian Progress Regression?

text: I know Gaussian Progress Regression is completed determined by its mean and covariance functions. For given training and test data, I can compute the covariance matrix correctly like a machine, and I can also describe the main procedure of GPR (compute the joint probability over training data and estimation point, cut a piece through the marginal likelihood....etc), but could anyone can describe what is Gaussion Progress for Regression using your OWN words?


title: What is an appropriate statistical test to identify significantly different timepoints in two timecourses?

text: I have two timecourses. Both are the same length. Both are univariate. Each represents the average EEG signal from a unique subgroup. The two subgroups do not have the same number of subjects. I would like to find individual timepoints where the difference between the two groups is statistically significant. It seems as though my options are:   * Perform a two-tailed t-test for each timepoint and correct for multiple comparisons.   * Perform an ANOVA analysis.   * Use a GLM to model each timecourse and then analyze the beta coefficients However, there seem to be drawbacks to each approach.   * Given the length of the EEG timecourse, and the relatively small difference between the two values, I worry about the risk of false negatives in the case of t-tests.   * ANOVA can tell me _if_ a given timepoint is different, but does not identify the actual timepoint of interest. Additionally, the assumption of independence would be violated because sequential measurements will be correlated.   * Linear analysis would not identify timepoints of interest. Any suggestions? Thanks for your time.


title: Conditional variance - $Var(X + U | X) = Var(U)$?

text: I am wondering if the following equality holds - $Var(X + U | X) = Var(U)$? where $X$ and $U$ are two independent random variables? It seems can we say $Var(X + U | X) = Var(X|X) + Var(U|X) = Var(X|X) + Var(U)$ as $U, X$ independent? And then $Var(X|X) = 0$? But that doesn't seem right. Can someone clarify this for me, I'm not very experienced with conditional variance.


title: Finding starting values for nls for critical exponential function

text: I'm trying to fit a critical exponential function using `nls`. The function is of the form $y = a+ (b + cx)r^x$ and has a single maximum/minimum, a single inflection point and an asymptote. I can sort of estimate starting values for the parameters from the data (eg the asymptote), and this works some of the times, but I'd like a more rigorous method to estimating the parameters prior to using `nls`. I found a suggestion here, which uses expected-value parameters to estimate the parameters for an exponential curve. I tried to do something similar, but the extra parameter makes things complicated. I also found a reparameterisation of the function which removes the `c` parameter by using expected-value parameter $$ y=[b(1-x/x_1)+(x/x_1)(y_1-a)/r^{x_1}]r^x+a $$ where `y1` corresponds to $x=x_1$, but I can't figure out how to estimate the remaining parameters. Any suggestions as to how to approach this would be appreciated! Some test data:               xvals=c(30.0, 46.0, 62.0, 78.0, 94.0, 110.0, 126.0, 142.0, 158.0, 174.0, 190.0, 206.0, 222.0, 238.0, 254.0)     yvals=c(-0.42286827, -0.6666351, -0.38212553, 0.10617201, 0.88309403, 0.97816225, 0.7471844, 0.3668751, 0.09183945, 0.11818037, 0.18587128, 0.22363502, 0.32338324, 0.27689159, 0.32192359)      Thanks!


title: Discounted cumulated gain

text: I've a little question regarding the Discounted Cumulated Gain (DCG) (Sorry, I couldn't find the papers of Järvelin and Kekäläinen). Can this evaluation- metric be used when a information retrieval system has only a binary classification of relevance? If this topic is misplaced, please tell where I should post it. Thank you


title: Total least squares curve fit problem

text: I am trying to fit a quadratic curve across a scatter plot of two variables. Since both variables are noisy I cannot use an ordinary least square regression (OLS) and I would like to have a symmetrical answer regardless which one I choose to be the independent variable. Total (orthogonal) least squares regression seem to be the answer. Unfortunately it seems that the two implementations that I have tried to use so far give me very bad fits. In the figure below the blue line is the OLS fit, which obviously could be improved. The red is the TLS fit using the matlab/octave code below which seems to be the standard approach using single-value decompositions (SVD). I copied it from the corresponding Wikipedia article here. I can see how it is trying to fit the upper-right denser cloud of points to the (very large) expense of the long tail. I have tested the code with simulated data and it works so is does not seem to be a bug. I have also tried this other matlab package that implements a different iterative approach but unless I give it very tight bounds to the model parameters also gives me "crazy" fits. I want to be able to perform TLS between many variables with many bootstrapping repeats to obtain confidence intervals and I cannot afford to come out with bounds for all those repeats eyeballing each graph; that would also tamper with the confidence interval estimation I guess. I reckon the problem is that the data is a bit too noisy but I can do little about it. You can find the data here. Any ideas how to go around this issue? Perhaps some know a more robust TLS method with an easy to use implementation available? There seem to be many academic articles but little runnable code. Many thanks in advance! V. ![enter image description here](http://i.stack.imgur.com/LFaO0.jpg)               function B = tls(xdata,ydata)     % xdata = nx3 matrix [x.^2 x ones(size(x))]     % ydata = nx1 matrix      m       = length(ydata);       %number of x,y data pairs     X       = [xdata];     Y       = ydata;     n       = size(X,2);          % n is the width of X (X is m by n)     Z       = [X Y];              % Z is X augmented with Y.     [U S V] = svd(Z,0);           % find the SVD of Z.     VXY     = V(1:n,1+n:end);     % Take the block of V consisting of the first n rows and the n+1 to last column     VYY     = V(1+n:end,1+n:end); % Take the bottom-right block of V.     B       = -VXY/VYY;     end


title: Two-way ANOVA repeated measures

text: I would like to know, if you got a significant interaction in a two-way repeated measures ANOVA in SPSS, how can you visualize the interaction (especially using syntax)?


title: Calculating the Similarity of Survey Responses

text: I was wondering if anyone had experimented with different functions for calculating the similarity of two sets of survey responses. I am going to be plugging it into a hierarchical clustering algo and need some sort of way to determine similarity. I was thinking of normalizing and using hyperbolic tangent. Or, coming up with my own function that essentially tells you how many shifts of the question responses it would take to arrive at the current set of responses. Any thoughts? Paul


title: Is there any reason to prefer a bayesian model with few variables?

text: I have two alternative hierachical bayesian models that were designed to the describe the same process (from a high-level point-of-view). Both model provides comparable (but not identical) inferences on test data-sets and both exhibits satisfactory behaviours. So from a practical point-of-view, I have no reason to prefer one model over the other (while they rely on different views). One important difference between them consists in the number of auxiliary latent rv composing the model: one model is composed $N+1$ (scalar) latent rv whereas the other is composed of (scalar) $2N+1$ latent rvs($N$ being the nber of observed sample). Is there any reason to prefer a model with less auxiliary variables (regarless of the computational time)?


title: What theory of estimation does apply when estimation follows model selection?

text: Almost two decades ago, Chatfield JRSS(1995)[vol.158,p.441] wrote that               Least-squares theory is known to not apply when the same data are used to formulate and      fit a model so that estimation follows model selection. Substantial model selection      biases can arise, particularly with subset selection methods in multiple regression      and in time series analysis. Unfortunately ways of overcoming the problem are not so clear.      At least for the ordinary regression model with normally distributed errors, the above problem applies to the maximum likelihood estimator. At least for certain Gaussian AR processes, conditional likelihood methods lead to least- squares problems. So the statement is quite general. Has any real progress been done in this direction over the last 20 years? If yes, what theory does apply, particularly in multiple regression and in time series analysis, when estimation follows model selection?


title: What’s wrong with this way of fitting time-dependent coefficients in a Cox regression?

text: I have a Cox proportional hazards model. Judging by Schoenfeld residual vs. time plots and corresponding tests for zero slope, there is clear violation of the PH assumption for several of the variables (too many to stratify). If we take the single-covariate case for simplicity, the original model looks like: $$ \hat\lambda(t) = \lambda_0(t)*\exp(X_1\hat\beta_{X_1}) $$ Given the violation of the PH assumption, I want to try fitting time-dependent coefficients for the problem variables. It seems intuitive to fit this: $$ \hat\lambda(t) = \lambda_0(t)*\exp(X_1\hat\beta_{X_1}+X_1t*\hat\beta_{X_1t}) $$ which simplifies to $$ \hat\lambda(t) = \lambda_0(t)*\exp(X_1(\hat\beta_{X_1}+t*\hat\beta_{X_1t})) $$ So, effectively, this would parametrically fit a coefficient $\hat\beta^\prime_{X_1}$ that is a linear function of time: $$ \hat\beta^\prime_{X_1} = \hat\beta_{X_1}+t*\hat\beta_{X_1t} $$ Before I could happily dance off into the sunset with this model, I learned it's invalid (Terry Therneau #1 and #2). But WHY is it invalid? Why do we have to construct a dataset with start-stop times instead?


title: System missing data

text: I'm currently struggling with my data. My dependent variable has an N of 332, out of which 180 are system missing, because of the routing in the questionnaire: First question: Have you taken in the past two years some substantial new machines into production? Second question: How many percent of your budget did they account for? So this means if people answered 'no' at the first question, they won't have answered the second question, and so they appear in the data list as system missing. Do I need to do some imputations, like 0%, or code it as -99, or just leave as is? What would you recommend? I want to test some hypotheses with this dependent variable through a multiple regression.


What does it mean to determine the power of your experiment via simulation? Does it mean to run the experiment over and over again and count how many times a specific set of null hypotheses are rejected given that they are false?


title: How to model multivariate input and multivariate output?

text: Given a set of multivariate inputs and corresponding multivariate output for each input, what statistical tools can be used to develop a model between multivariate input and output? I am wondering if using probabilistic graphical models is the way to go. I have the detailed values of input and output simultaneously and such 50 input-output pair. Both input and output are around 100-dimensional vector (I know this still sounds vague, but I not permitted to discuss the details of input and output sequences because of confidentiality issues). The question I am trying to answer is "Given the input vector, can I develop a probabilistic model to determine the output sequence?." This problem is relatively new to me, so initially I had formulated it wrongly by stating that the inputs and outputs are histograms. I have edited it now to avoid any confusion.


title: Hypothesis Testing with Spatial Data

text: I have production data for several wells. I know that a subset has been treated with a product, and I'd like to determine if the difference in production is due to the product treatment, or the placement of the wells. For each test well I've collected the $K$ nearest neighbors required to obtain $N$ control wells. This gives me $N$ control wells, and $K-N$ test wells. I then perform a two-tailed z-test, using the standard deviation of the production of the $K$ wells as the sample standard deviation. I was wondering if this sounds like a decent approach to this problem, or if there is a more sophisticated way of incorporating spatial data for testing the hypothesis that the production for the test and control groups is the same.


title: hierarchical logistic regression Block non-significant, interaction significant

text: I'm using hierarchical logistic regression and having some difficulty. Im trying to predict self-harm vs none based on IVs age, gender, alcohol, educational attainment. They are all yes/no except educational attainment (e.g. low, mid, high). I have added age, gender, educational attainment in step 1, and alcohol in step 2, step three is my interaction alcohol*educational attainment. all are significant except educational attainment. While the model is significant for model 3, the Block is ns (p=.091). I think this is because while alcohol*edu attainment (low vs high) is significant (OR 6.655, p=.038), alchol*edu attainment (low vs medium) is highly non-significant (OR=.179, p=.282). The issue is that i believe the interaction is really important. When i run the model again without the 'middle' educational attainment the Block becomes significant. Is the reason behind this believe it to be? Is there anyway i can justify using the model with the interaction despite the fact Block is non-significant?


title: Problem simulating AR and MA models using filters

text: I do not know how to use filter to simulate AR and MA models. To me it looks the same way for both MA and AR Estimate AR so then how do I know that the model is AR or MA ? For example, **Problem1:** For the MA(4) model the constant term is 0.5 and the coefficients are `$h_{MA} = [0.8 0.0 0.0 0.2]$`. The innovation process $\epsilon$ is a deterministic signal Logistic Map. The functional form of the MA(4) model looks like: $y(t)= 0.5 + \epsilon(t)+ 0.8\epsilon(t-1) + 0.2\epsilon(t-4) + additive white Gaussian noise$ This is how I generated MA(4) model. I am having difficulty in including the $\epsilon(t)$ in the filter.                x(1) = 0.1;      N = 100;         for j =1 : N                   x(j+1) = 4*x(j)*(1-x(j));         end          y1 = filter(h,1,x);     noise = 10; %SNR in Db     y = awgn(y1,noise,'measured');      Is the method correct? **Problem2:** How do I generate an AR model based on the same logistic map and then add noise for an AR(p) for example $y(t) = a1y(t-1) + a2y(t-2) + \epsilon(t) + additive white Gaussian noise$? **Problem3:** When calculating error between desired coefficients and estimated, should the Mean Square error be (A) $sum((\hat{H}-h).^2)/5$ or (B) $sum((\hat{H}-h).^2)/4$ ? i.e should I consider the constant term or not?


title: Why is this a good approximation of cross entropy?

text: Given a random variable $X=x_i$ for $i=1,...,n$, with true distribution $p(x)$ and approximate distribution $q(x)$, its cross entropy is given by $$H(p,q) = -\sum_{i=1}^np(x_i)\log q(x_i).$$ However in practice we find that the true distribution $p(x)$ is rarely if ever known, and thus instead we compute an "approximation" to the cross entropy, given by $$-\sum_{i=1}^n\frac{1}{n}\log q(x_i).$$ Why is this a good approximation of the cross entropy? I can't figure out why it would be, even if we sample enough points from the true distribution to make $q(x)$ a really good approximation for $p(x)$, why would we want all weights to remain the same? In that case wouldn't we want to just replace $\frac{1}{n}$ with $q(x_i)$? And if $q(x)$ is not known to be a good approximation, then setting $p(x_i)=\frac{1}{n}$ seems like a completely arbitrary choice.


title: what is meant by binary RBM?

text: I understand what is RBM. I also accept that in binary RBM, states of RBM are binary. Is that also means the input datasets should be binary?


title: Communicating classes in a Markov Chain

text: If my state space is $E=1,2,3,4,5,6,7$ and transition matrix is:   \begin{pmatrix} .5 &0 &0 &.5 &0 &0 &0 \\ 0 &.2 &.3 &.1 &.4 &0 &0 \\ .4 &0 &0 &.6 &0 &0 &0 \\ 0 &0 &.2 &.8 &0 &0 &0 \\ 0 &.3 &0 &0 &0 &.3 &.4 \\ 0 &0 &0 &0 &.2 &.6 &.2 \\ 0 &0 &0 &0 &1 &0 &0 \end{pmatrix} can someone tell me the communicating classes, closed sets and classify all states? I think the communicating classes are $\{1,3,4\},\ \{2\},\ \{5,6,6\}$, and that the closed set is $\{1,3,4\}$.


title: Sample probability for binomial distribution

text: $X\sim B(40, 0.3), \bar X$ is the sample mean, and $n$ is number of observations. I need help in finding $P(\bar X\ge13)$ for a sample of $49$. I have tried the following, please guide. $\mu=40 \times 0.3=12,$ and $\sigma =\sqrt{40\times0.3\times0.7}=2.8982$ by standardizing I got $P(Z>\dfrac{13-12}{\frac{2.8982}{7}})=1-\phi(2.4152922)=0.0078612~$ correct answer should be 0.0084. Please help Also I am having difficulties in finding the smallest value of $n$ for which $P(\bar X\ge13)\lt0.001$, where $n$ is the number of observations. Any help is highly appreciated.


title: How does one derive the survival function from an accelerated time failure model in R?

text: I'm currently working with accelerated time failure models. However, when I fit a model using `survreg` in R, the `predict` function only seems to return expected lifetime. I'm wondering how one can use a fitted accelerated time failure model to predict probability of survival within a given time interval, conditional on covariates and time already survived. That is, I would like to know: ![enter image description here](http://i.stack.imgur.com/mO35W.jpg) Ideally, I would have a predict function that accepts a fitted model object, the times in question, and new data containing covariates as arguments. The function then returns probability of survival up to each time, i.e. the conditional survival function, for each observation in the new dataset. **Edit** To make this reproducible:               lfit <- survreg(Surv(time, status) ~ ph.ecog, data=lung)     pct <- 1:98/100 # The 100th percentile of predicted survival is at +infinity     ptime <- predict(lfit, newdata=data.frame(ph.ecog=2), type='quantile',                  p=pct)      `ptime` contains the number of days associated with each percentile. I want the inverse: to supply a time value and get the probability of survival.


title: How often to update a forecasting model

text: If I have a data set with daily, weekly and annually seasonality, how often should I update my forecasting model? As I have heard, forecasting models out there can have a good prediction for up to a day ahead for these types of data, hence should I update my model daily to find new values for seasonality for the upcoming periods? Thanks


title: Hypothesis test for whether correlation coefficient is greater than specified value?

text: I want to test whether a correlation coefficient is greater than 0.5.   * What should the null and alternative hypothesis be?   * Should I do a two tailed test or one tailed test? Thanks